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Abstract 

The objective of this research is to improve and implement the filtered mass 
density function (FDF) methodology for large eddy simulation (LES) of high- 
speed reacting turbulent flows. We have just completed three years of this 
research. This is the Final Report on our activities during the period: Septem- 
ber 1, 1999 through August 31, 2002. 

Dr. J. Philip Drummond (Hypersonic Propulsion Branch, NASA LaRC, 

Mail Stop 197, Tel: 757-864-2298) is the Technical Monitor of this Grant. 


1 Introduction 

An issue of current interest to NASA is associated with design of various components 
involved in air-breathing propulsion systems such as the scramjet [1]. A successful 
design needs the resolution of numerous issues such as fuel injector design and opti- 
mization, flame holding associated with the injection, critical design studies of certain 
areas (not the whole integrated engine), etc. Resolving these issues is of crucial im- 
portance to achieve the final goal of producing a more optimum aerospace propulsion 
system for hypersonic vehicles. To achieve this goal, however, there is a demand for 
development of robust tools that can aid in the design procedure. This is due to the 
fact that the physics of high-speed reactive flows is rich with many complexities. Few 
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examples of the physical issues of current interest are the questions associated with 
the chemical and thermodynamical non-equilibrium effects, the cause and effect of 
turbulence, the interaction of turbulence and chemistry, the real gas effects at high 
temperatures, etc. 

The need for advanced computational methods for the analysis of high speed propul- 
sion systems is obvious [1]. Within the past three decades the NASA Langley Research 
Center (LaRC) has been at the forefront in making use of advanced computational 
methods for investigations of such systems. Large eddy simulation (LES) is regarded 
as one the most promising means of simulating turbulent reacting flows [2,3]. Amongst 
the various LES strategies, the approach based on the probability density function 
(PDF) has proven particularly effective [4-19]. The formal means of conducting such 
LES is by consideration of the “filtered density function” (FDF) which is essentially 
the filtered fine-grained PDF of the transport quantities [20]. Colucci et al. [8] devel- 
oped a transport equation for the FDF in constant density turbulent reacting flows. 
Jaberi et al. [9] extended the methodology for LES of variable density flows by con- 
sideration of the “filtered mass density function” (FMDF) which is the mass weighted 
FDF. The fundamental property of the PDF methods is exhibited by the closed form 
nature of the chemical source term appearing in the transport equation governing the 
FDF (FMDF). This property is very important as evidenced in several applications 
of FDF for LES of a variety of turbulent reacting flows [8-12,16]. 


2 Accomplishments 

The goal of this research was to improve the capabilities of the FDF method and to 
implement it for LES of chemically reacting turbulent flows. We feel that we have 
been very successful in achieving the specific objectives of this work. These objectives 
were: 

1. Development and implementation of the velocity filtered density function 
(VFDF) for LES. 

2. Development and implementation of the joint velocity-scalar filtered den- 
sity function (VSFDF) for LES. 

3. Implementation of the FDF for LES of hydrocarbon diffusion flames and 
comparison with experimental data. 
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The efforts pertaining to the first two objectives are completed and are fully described 
in journal articles (included here as Appendix I and Appendix II). The work pertain- 
ing to the third objective is premature for publication. Also, during the work on 
this project, the PI was invited to deliver a Keynote Lecture at the Third AFOSR 
International Conference on DNS and LES (Arlington, TX, August 5-9, 2001). A 
copy of the review article on this tutorial lecture is given in Appendix III. The efforts 
pertaining to each of these objectives are summarized in the next three subsections. 


2.1 Velocity Filtered Density Function for Large Eddy Sim- 
ulation of Turbulent Flows 

In this part of our work, a methodology termed the “velocity filtered density function 
(VFDF) is developed and implemented for large eddy simulation (LES) of turbulent 
flows. In this methodology, the effects of the unresolved subgrid scales (SGS) are 
taken into account by considering the joint probability density function (PDF) of all 
of the components of the velocity vector. An exact transport equation is derived for 
the VFDF in which the effects of the SGS convection appear in closed form. The 
unclosed terms in this transport equation are modeled. A system of stochastic differ- 
ential equations (SDEs) which yields statistically equivalent results to the modeled 
VFDF transport equation is constructed. These SDEs are solved numerically by a 
Lagrangian Monte Carlo procedure in which the Ito-Gikhman character of the SDEs 
is preserved. The consistency of the proposed SDEs and the convergence of the Monte 
Carlo solution are assessed by comparison with results obtained by an Eulerian LES 
procedure in which the corresponding transport equations for the first two SGS mo- 
ments are solved. The VFDF results are compared with those obtained via several 
existing SGS closures. These results are also analyzed via a priori and a posteri- 
ori comparisons with results obtained by direct numerical simulation (DNS) of an 
incompressible, three-dimensional (3D), temporally developing mixing layer. 

This work is fully described in Appendix II. 
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2.2 Velocity-Scalar Filtered Density Function for Large Eddy 
Simulation of Turbulent Flows 

In this part of our work, a methodology termed the “velocity-scalar filtered density 
function” (VSFDF) is developed and implemented for large eddy simulation (LES) 
of turbulent flows. In this methodology, the effects of the unresolved subgrid scales 
(SGS) are taken into account by considering the joint probability density function 
(PDF) of the velocity-scalar field. An exact transport equation is derived for the 
VSFDF in which the effects of the SGS convection and chemical reaction are closed. 
The unclosed terms in this equation are modeled. A system of stochastic differ- 
ential equations (SDEs) which yields statistically equivalent results to the modeled 
VSFDF transport equation is constructed. These SDEs are solved numerically by a 
Lagrangian Monte Carlo procedure in which the Ito-Gikhman character of the SDEs 
is preserved. The consistency of the proposed SDEs and the convergence of the Monte 
Carlo solution are assessed by comparison with results obtained by a finite difference 
LES procedure in which the corresponding transport equations for the first two SGS 
moments are solved. The VSFDF results are compared with those obtained via other 
SGS closures, and all the results are assessed via comparison with data obtained by 
direct numerical simulation (DNS) of a temporally developing mixing layer involv- 
ing transport of a passive scalar. It is shown that the values of both the SGS and 
the resolved components of all second order moments including the scalar fluxes are 
predicted well by VSFDF. The sensitivity of the model’s (empirical) constants are 
assessed and it is shown that the magnitudes of these constants are in the same range 
as that typically employed in PDF methods. 

This work is fully described in Appendix II. 


2.3 Scalar Filtered Density Function for Large Eddy Simu- 
lation of a Diffusion Flame 

We have started work on the use of the filtered density function (FDF) method for 
large eddy simulation (LES) of the piloted jet flame configuration as considered in 
the experiments of the Combustion Research Facility at the Sandia National Lab- 
oratory [21-25]. This flame has been the subject of broad investigations by other 
computational/modeling methodologies [24], In the experiments, three basic flames 
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are considered, identified by Flames D, E, and F. The geometrical configuration in 
these flames is the same, but the jet inlet velocity is varied. In Flame D, the fuel jet 
velocity is the lowest and the flame is close to equilibrium. The jet velocity increases 
from flames D to E to F, with noticeable non-equilibrium effects in the latter two. 

Some preliminary work has already been completed in which the LES/SFDF is used 
for predictions for the Sandia piloted jet flame. In these simulations, combustion is 
modeled via two chemistry models: (1) an equilibrium model via realistic kinetics, 
(2) a finite rate, single-step model for non-equilibrium flames. In (1), the LES/SFDF 
is employed in conjunction with equilibrium methane-oxidation model. This model 
is enacted via “flamelet” simulations which consider a laminar counterflow (opposed 
jet) flame configuration [26-29]. The full methane oxidation mechanism of the Gas 
Research Institute (GRI) [30,31] accounting for 53 species and 325 elementary re- 
actions is employed. At low strain rates, the flame is close to equilibrium. Thus, 
the thermo-chemical variables are determined completely by the “mixture fraction.” 
This flamelet library is coupled with our LES/SFDF solver in which transport of the 
mixture fraction is considered. It is useful to emphasize here that the PDF of the mix- 
ture fraction is NOT “assumed” a priori (as done in almost all other flamelet based 
LES [32—42] ) ; rather, it is calculated explicitly via the SFDF. In (2), the methane 
oxidation is modeled via a finite-rate, single-step kinetics model [43], The system 
of nonlinear ordinary different equations representing reactant conversion is solved 
via LES/SFDF for all of the scalar variables (mass fractions and enthalpy). With 
this, some aspects of non-equilibrium chemistry are taken into account, albeit in an 
idealized manner. 

With the equilibrium chemistry model, the results obtained by LES/SFDF/flamelet 
are compared with Sandia’s experimental data [21-25]. Up to now, we have only 
simulated Flame D and we have observed good agreements. But more work is required 
to ensure the accuracy of our results. Simulations of flames E and F have not been 
yet attempted. 

In closure of this section, we would like to state that since our early work on FDF, 
this methodology has been used by several other investigators (e.g. Refs. [12,16,19]). 
Please see Appendix III for a recent review of current state of progress in LES/FDF. 
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3 Enhancement of Technology and Education 


With completion of this research, we feel that we have been able to contribute to 
maintain U.S. leadership in a technology which is ot significant importance to NASA. 
This is a very important matter, particularly when considering the extent of ongoing 
concentrated efforts in Europe and Asia in the efforts proposed here. To appreci- 
ate the seriousness of the competition, we indicate that the number of European 
researchers who are currently involved in mathematical & computational modeling 
of turbulent reacting flows is significantly larger than that in the U.S. In fact, con- 
sidering only England, France and Germany, the number of senior investigators (not 
including graduate students) who are working on the specific research field of “LES of 
turbulent combustion” in these three countries is larger than that in the whole U.S.! 
This indicates the wide recognition of the importance of this research field in other 
countries. We feel that the approach proposed here is very promising in dealing with 
this challenging problem in addition to enhancing the current state of knowledge in 
turbulent combustion. 

In order to demonstrate our visibility in this research, here we shall list all the awards 
and some of noticeable achievements of the personnel involved in this program. 


3.1 Graduate Students 

Involvement of students in research is an issue which is taken very seriously at our 
University We are committed to recruiting excellent quality students and involving 
them in high caliber research. During the past three years, the following students 
have been supported under this Grant. 

1. Dr. Laurent Y.M. Gicquel (Ph.D. 2001). Currently: Research Scientist, CERFACS, 
Toulouse, France. NASA’s support is acknowledged in Ph.D. Disseration [44]. 

2. Mr. Tomasz G. Drozda (M.S. 2002). Currently: Ph.D. candidate at the University 
of Pittsburgh. NASA’s support is acknowledged in M.S. Thesis [45] 

3. Mr. M. Reza H. Sheikhi (M.S. 2002). Currently: Ph.D. candidate at the University 
of Pittsburgh. 
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3.2 Awards and Honors 


1. Peyman Givi: Promoted to University at Buffalo Distinguished Professor (2002). 

2. Peyman Givi: Received Professor of the Year Award, Tau Beta Pi Engineering 
Honor Society, New York Nu Chapter (2001-2002). 

3. Tomasz G. Drozda: Second Prize winner at the Graduate Technical Paper Com- 
petition of AIAA Northeast Regional Student Conference, Rensselaer Polytechnic 
Institute, Troy, NY. Title of paper: “A Hybrid Stochastic-Deterministic Methodol- 
ogy for Large Eddy Simulation of Scalar Mixing and Reaction in Turbulent Flows,” 
April, 2002. 

4. Profile Featured in: 

• Reporter: 13 named UB Distinguished Professor 33 (28), p. 1, May 9, 2002. 

• The Buffalo News: Rolls-Royce using UB technique to refine its engines, p. E4, 
March 19, 2002. 

• SEAS News: LES results equal more expensive supercomputer simulations, 

VIII(l), p. 6, Spring (2002). 

• Reporter: New method produces “super” results. 33(21), p. 1, March 14, 2002. 

• Reporter. Graduates bringing recognition to lab, 32(20), p. 4, February 15, 

2001 . 

• ASEE Recruitment Video , the NASA Langley Research Center, Office of Edu- 
cation, Hampton, VA, 1999-2000. 


3.3 Publications 

In all of the following publications, the support from NASA Langley is greatfully 
acknowledged: 

Tnvited Review Articles: 

1. P. Givi, “A Review of Modern Developments in Large Eddy Simulation of Tur- 
bulent Reactive Flows,” Chapter in DNS/LES: Progress and Challenges , pp. 81-92, 
Editors: C. Liu, L. Sakell and T. Beutner, Greyden Press, Columbus, OH, 2001. 
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2. F.A. Jaberi, F. Mashayek, C.K. Madnia, D.B. Taulbee and P. Givi, “Advances 
in Analytical Description of Turbulent Reacting Flows,” Chapter 9 in Advances in 
Chemical Propulsion , pp. 149-164, Editor: Gabriel, D. Roy, CRC Press LLC, Boca 
Raton, FL, 2002. 

Refereed Papers: 

1. M.R.H. Sheikhi, T.G. Drozda, P. Givi, and S.B. Pope, “Velocity-Scalar Fil- 
tered Density Function for Large Eddy Simulation of Turbulent Flows,” submitted 
to Physics of Fluids (2002). 

2. L.Y.M. Gicquel, P. Givi, F.A. Jaberi, and S.B. Pope, “Velocity Filtered Density 
Function for Large Eddy Simulation of Turbulent Flows,” Physics of Fluids , 14(3), 
1196-1213, 2002. 

3. L.Y.M. Gicquel, P. Givi, F.A. Jaberi, and S.B. Pope, “Velocity Filtered Density 
Function for Large Eddy Simulation of a Turbulent Mixing Layer, in DNS/LES: 
Progress and Challenges , pp. 327-334, Editors: C. Liu, L. Sakell and T. Beutner, 
Greyden Press, Columbus, OH, 2001. 

Conference Papers: 

1. M.R.H. Sheikhi, T.G. Drozda, P. Givi, and S.B. Pope, “Velocity-Scalar Filtered 
Density Function for Large Eddy Simulation of Turbulent Flows, submitted for 
presentation at the 55th Annual Meeting of the Division of Fluid Dynamics of the 
American Physical Society, Dallas, TX, November 2002. 

2. P. Givi, L.Y.M. Gicquel, F.A. Jaberi and S.B. Pope, “PDF Methods for Large 
Eddy Simulation of Turbulent Reactive Flows,” Proceedings of IUTAM Symposium 
on Turbulent Mixing and Combustion, pp. 96-97, Kingston, Ontario, Canada, June 
2 - 6 , 2001 . 

3. L.Y.M. Gicquel, P. Givi, F.A. Jaberi and S.B. Pope, “Velocity Filtered Density 
Function for Large Eddy Simulation of Turbulent Flows,” Bulletin of the American 
Physical Society , 45(9), p. 129, 53rd Annual Meeting of the Division of Fluid Dy- 
namics of the American Physical Society, Washington, DC, Nov. 19-21, 2000. 

4. F.A. Jaberi, S. James and P. Givi, “Large Eddy Simulations of Methane Jet 
Flames,” Bulletin of the American Physical Society , 45(9), p. 72, 53rd Annual Meet- 
ing of the Division of Fluid Dynamics of the American Physical Society, Washington, 
DC, Nov. 19-21, 2000. 
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5. S.C. Garrick, F.A. Jaberi and P. Givi, “Large Eddy Simulation of Turbulent 
Reacting Round Jets,” Bulletin of the American Physical Society , 44(8), p. 44, 5'2nd 
Annual Meeting of the Division of Fluid Dynamics of the American Physical Society, 
New Orleans, LA, November 21-23, 1999. 
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Abstract 

A methodology termed the “velocity filtered density function” (VFDF) is developed 
and implemented for large eddy simulation (LES) of turbulent flows. In this method- 
ology, the effects of the unresolved subgrid scales (SGS) are taken into account by 
considering the joint probability density function (PDF) of all of the components of 
the velocity vector. An exact transport equation is derived for the VFDF in which 
the effects of the SGS convection appear in closed form. The unclosed terms in this 
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transport equation are modeled. A system of stochastic differential equations (SDEs) 
which yields statistically equivalent results to the modeled VFDF transport equation 
is constructed. These SDEs are solved numerically by a Lagrangian Monte Carlo pro- 
cedure in which the Itd-Gikhman character of the SDEs is preserved. The consistency 
of the proposed SDEs and the convergence of the Monte Carlo solution are assessed 
by comparison with results obtained by an Eulerian LES procedure in which the cor- 
responding transport equations for the first two SGS moments are solved. The VFDF 
results are compared with those obtained via several existing SGS closures. These re- 
sults are also analyzed via a priori and a posteriori comparisons with results obtained 
by direct numerical simulation (DNS) of an incompressible, three-dimensional (3D), 
temporally developing mixing layer. 

Corresponding Author: Peyman Givi, Tel: (716) 645-2593 (x. 2320). Fax: 716-645- 
3875. E-mail: givieng.buffalo.edu. 


2 



1 INTRODUCTION 


The probability density function (PDF) approach has proven useful for large eddy 
simulation (LES) of turbulent reacting flows . 1-15 The formal means of conducting such 
LES is by consideration of the “filtered density function” (FDF) which is essentially 
the filtered fine-grained PDF of the transport quantities. In all previous contribu- 
tions, the FDF of the “scalar” quantities is considered: Gao and O’Brien , 3 Colucci et 
al ., 6 Reveillon and Vervisch 7 and Zhou and Pereira 12 developed a transport equation 
for the FDF in constant density turbulent reacting flows. Jaberi et al . 8 extended 
the methodology for LES of variable density flows by consideration of the “filtered 
mass density function” (FMDF), which is essentially the mass weighted FDF. The 
fundamental property of the PDF methods is exhibited by the closed form nature 
of the chemical source term appearing in the transport equation governing the FDF 
(FMDF). This property is very important as evidenced in several applications of FDF 
for LES of a variety of turbulent reacting flows . 6-10,12 However, since the FDF of only 
the scalar quantities are considered, all of the “hydrodynamic” effects are modeled. In 
all previous LES/FDF simulations, these effects have been modeled via “non-FDF” 
methods. 

The objective of the present work is to extend the FDF methodology to also include 
the subgrid scale (SGS) velocity vector. This is facilitated by consideration of the joint 
“velocity filtered density function” (VFDF). With the definition of the VFDF, the 
mathematical framework for its implementation in LES is established. A transport 
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equation is developed for the VFDF in which the effects of SGS convection are shown 
to appear in closed form. The unclosed terms in this equation are modeled in a fashion 
similar to those in the Reynolds-averaged simulation (RAS) procedures. A Lagrangian 
Monte Carlo procedure is developed and implemented for numerical simulation of the 
modeled VFDF transport equation. The consistency of this procedure is assessed by 
comparing the first two moments of the VFDF with those obtained by the Eulerian 
finite difference solutions of the same moments transport equations. The results of the 
VFDF simulations are compared with those predicted by the Smagorinsky 16 closure, 
and the “dynamic” Smagorinsky model. 17-19 The VFDF results are also assessed via 
comparisons with direct numerical simulation (DNS) data of a three-dimensional (3D) 
temporally developing mixing layer in a context similar to that of Vreman et al. 20 

This work deals with LES of the velocity field in a constant density, non-reacting 
flow. Consideration of the joint velocity-scalar FDF (or FMDF) in variable density, 
chemically reacting flows will be the subject of future work. It is in this context that 
the approach has its principal advantage: convective transport (of momentum and 
species) is in closed form. 


2 FORMULATION 

In the mathematical description of incompressible (unit density) turbulent flows, the 
primary transport variables are the velocity vector, «i(x,t) (z = 1,2,3), and the 
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pressure, p(x,t), field. The equations which govern transport of these variables in 
space (xj) and time ( t ) are 


du z 

dxi 


- 0 , 


duj duiUj dp da, 

dt dxi dxj dxi 


ij 


For a Newtonian fluid, the viscous stress tensor is represented by 


'1J 


= V 



( 1 ) 


( 2 ) 


where v is the kinematic viscosity and is assumed constant. 

Large eddy simulation involves the spatial filtering operation 21 23 

/ + oo 

/(x',t)£(x',x)dx\ (3) 

-00 

where Q denotes the filter function, ( f(x,t)) L represents the filtered value of the 
transport variable /(x, t), and f = f — {/)l denotes the fluctuations of / from the 
filtered value. We consider spatially and temporally invariant and localized filter 
functions, thus £(x',x) = G(x' - x) with the properties, 21 G(x) = G(-x), and 
f^° oo G(x)dx = 1. Moreover, we only consider “positive” filter functions 24 for which 
all the moments f^ x x m G(x)dx exist for m > 0. The application of the filtering 
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operation to the instantaneous transport equations yields 


d{uj) L _ Q 

dxi ^ 

d(uj) L = d(p) L d{a l3 ) L _ dT L (ui,Uj) 

dt dxi dxj dxi dxi 

where Ti(ui,u 3 ) = ( UiU 3 )i — (u{) l(u 3 ) l denotes the “generalized SGS stresses . 18 
These stresses satisfy 18 


d 0 

[TL(Ui, U 3 )] + [(Uk) L t l( u u u j)] = — 


dT l3 k 

dx k 


— n 0 - + Pa — £ 


u 


-y 


( 5 ) 


In this equation, T i3 , t = Ti(ui,Uj,Uk) — is the SGS turbulent trans- 

port tensor where T L (ui,Uj,Uk) = (uiUjU^L — ( Ui)LT~L(uj,Uk ) — {uj) iTi(ui,Uk) — 
{u k )LT L {ui,u 3 ) - (u t )L(^j>L<UA')L - 18 The other terms are the SGS pressure-velocity 
scrambling tensor, §£:), the SGS production rate tensor, 

p i3 = - T L (tn,u k and the SGS dissipation rate tensor, 
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3 VELOCITY FILTERED DENSITY FUNCTION (VFDF) 


A Definitions 


The “velocity filtered density function” (VFDF), denoted by Pl , is formally defined 


/ +OO 

Q [v» w(x', t)] G{x! - x)dx', 

-OO 

3 

Q [v, u(x, t )] = 5\v - u(x, t)} = - Ui(x, <)] 


where 5 denotes the delta function and v is the velocity state vector. The term 
g[v,u(x,t)] is the “fine-grained” density, 25,26,11 and Eq. (6) defines the VFDF as the 
spatially filtered value of the fine-grained density. With the condition of a positive 
filter kernel, 24 P L has all the properties of the PDF. 26 For further developments, it is 
useful to define the “conditional filtered value” of the variable Q(x, t) by 


(Q{x,t)\u(x,t) = v) L = 0 Q\v) l = 


” QiA *)g[«, m(x', *)] Gjx^_ - x)dx ; , 

Pl( v 'i x , t) 
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where ( a\(3)i denotes the filtered value of a conditioned on (3. Equation (7) implies 


(i) For Q{x, t) = c, (Q(x,t)\v) L = c, 

{ii) For Q(x, t) = Q(u(x,t)), (Q(x, t)\v) L = Q{v), (8) 

/ +oo 

(Q(x, 0 I v )l-Pl(v; x, t)dv, 

■OO 

where c is a constant, and Q(x, t) = Q(u(x, t )) denotes the case where the variable Q 
is completely described by the variable u(x, t). From these properties it follows that 
the filtered value of any function of the velocity variable is obtained by integration 
over the velocity space 


(<2(x,t)) L 



Q(v)P L (v,x,t)dv. 


( 9 ) 


B VFDF Transport Equation 


The exact transport equation for the VFDF is derived in this section. Two forms of 
this equation are considered similar to those previously developed in PDF methods. 27 ' 31 
The starting point is to consider the time-derivative of Eq. (6), 


gP L (u;x,t) 

dt 



dg [u,u(x ; ,£)] 
dt dvi 


G(x' - 



dvj^t) 

dt 


p[v,u(x',t)] G(x' 


x) dx' 

— x) dx' . 


( 10 ) 
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This combined with Eq. (7) yields 


dP L (v,x.,t) 

dt 


d / dui 

d^ t [\ ~dt V 


P L (v,x,t) . 

L 


(ii) 


Substituting Eq. (1) into Eq. (11) yields, 


dPx,(p;x,t) _ d 

dt dvi 




+ 

L 



\v 


P L {v,x,t) 


L J 


(12) 


With the relation 


d_ 

dx Vi 


dujU k i 
dx k 


v 


P L (v;x,t) 

L 


= ~Vi fc 


ap L (^;x,o 

dx k 


(13) 


and decompositions 


v k P L = {u^lPl + [vk - 


' d P I, A p _ 9 (p)l p 

q — |v ) r L — — -Pl + 

^ OX{ f l C/Xj 


dp 


W 


9 (p)i 


dxi'~ / L dxi 


Pl , 


/da ikl \ d(a lk ) L \/da ik ^\ d(a ik ) L 

Wv/ L = —x, r Pl + *rP) L ’ “*r 


Pl , 


(14) 


the VFDF transport equation becomes, 


DPl 

Dt 


+ 


dx 


d r / , np , d(p) L dP L d(a lk ) L dP L 

[(v k — {U k )ijPi\ + 


dxi dv. 


±\((^-\v 

dvi \ \ dxi 


diph 

dx. 


--[( 

Sui [\ 


dx k dvi 
' /da. 


ik i 


d(a lk )i 


\dx k ^ V / L dx k ^ Pl 


, (15) 
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where ^ ^ + {u k ) L g denotes the “filtered” material derivative. 

Equation (15) is an exact transport equation for the VFDF. The first term on the 
right hand side represents the SGS convection of the VFDF in physical space and 
is closed. The second and third terms (which are also in closed form) represent 
the convection in velocity space due to the resolved pressure gradient and molecular 
diffusion, respectively. The last two terms are unclosed and denote convective effects 
in the velocity space due to SGS pressure gradient and SGS diffusion. 

Alternatively, the conditional diffusion term in Eq. (15) can be represented as: 


d_ 

dvi 



P L (u;x,t) 

L 


^ g 2 PL(^;x,t) 
dx k dx k 


d 2 

dvidvj 


dui dv,j 

V — 

dx k dx k 


\ v 


P L (v\x,t) , 
(16) 


in which the second term on the RHS involves the conditional expected dissipation. 
With this, the alternate form of the VFDF transport equation is: 


DPl 

Dt 


dx k 

d 


9 [(*-<» t ) L m+ mLdPL 


+ 


dv; 


dp 
dx i 


v ) - 


dx 

d_ {ph 

dxi 


dvi 


+ v 


d 2 P L 

dx k dx k 


d 2 

dvidvj 



dui duj 
dx k dx k 


|n 


Pl ■ 


(17) 


Equation (17) is another exact transport equation for the VFDF. The first term on 
the right hand side represents the SGS convection of the VFDF in physical space 
and is closed. The second term corresponds to the convection in the velocity space 
due to the resolved pressure gradient. The third term represents molecular diffusion 
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of the VFDF in physical space. The closure problem is associated with the last two 
terms. These represent, respectively, convection in velocity space by the unresolved 
SGS pressure gradient, and diffusion in velocity space by SGS dissipation rate tensor. 


C Modeled VFDF Transport Equations 

The generalized Langevin model (GLM) 27,32 is employed for closure of the VFDF 
transport equation. Here we introduce two modeled VFDF equations, which are 
denoted by “VFDF1” and “VFDF2.” These are presented in order. To close Eq. (17), 
VFDF1 is 


d_ ' 

dvi 




d 2 f/ duj_duj \ p 
dvidvj \ dx k dx k V / L L 


d 

dv t 


[ Gij (vj - ( Uj) L ) Pl} + -^C 0 £ 
+ &Pl 

dxk dxk dv^Vj 


d 2 P L 


dvidvi 


T 2v 


dfah 

dx k 


d 2 P L 

dx k dv x 


(18) 


To close Eq. (15), VFDF2 is: 


d_ ' 

dvi 



djph 

dxi 


d(Ji k 

dx k 


V 



d_ 

dvt 


Gij ( Vj - (uj) L ) Pl ] + r Cq £ 


&Pl 

dv t dvi 


(19) 


Note that these models (i.e., the first two terms on the right hand sides of Eqs. (18) 
and (19) are the same, but that they model slightly different quantities. With this 
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closure, the two terms in Gq and e jointly represent the SGS pressure-strain and SGS 
dissipation. These are modeled as: 26,11 


Gij = - Q + jC„) S,„ e = C, fc 3/ 7A t> ui = e/k. 


( 20 ) 


where uj is the SGS mixing frequency, is the filter width, k — is the 

SGS kinetic energy, and £ - \eh is the SGS dissipation rate. 

With the GLM, the two forms of the VFDF transport equation are, 


DPl 

Dt 


d 


{Vk 


+ 2 v 


dx k 
d{u t ) L d 2 P L 


(^*)z.) Pl 

d 


3 ^ + , 

I r\ r\ ' 


d 2 P L 


dxi dvi dx k dx k 


+ is 




dx k dx k dvi dv l 


. Gy j ( v j ~ ( u j)L ) Pl } + ^ C 0 e 


dx k 

&Pl 

dvidvi 


dx k 


&Pl 

dv t dvj 


(21) 


for VFDF1, and 


DPl 

Dt 


--S- [K - (' u k ) L )P L ] + 

dx k 


d&hdPL 

dxi dvi 


djoikh dPi 

dx k dvi 


d [ Gii 


Ml) Pl] + 2 Co £ 


&Pl 

dvidvi 


( 22 ) 


for VFDF2. Hereinafter, Eqs. (21) and (22) are referred to as “VFDF1” and “VFDF2,” 
respectively. The difference between these two equations is in the different treatment 
of the closed viscous terms. 
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D Transport Equations for Moments 


The zeroth, first and second moment equations corresponding to these two formula- 
tions are: 


For VFDF1: 


dfah _ n 

dxi 

d^Ujh d(ui) L {uj) L _ d(p) L + v d 2 {u 3 ) L _ dT L ( Ul ,u 3 ) 

dt dx l dx 3 dxidxi dx t 


(23) 


~[T L (Ui,U])} + U k ) L ri{Ui,Uj )] = -J^-[TL{Ui,Uj,U k ) - V-^[T L {Ui,Ujj\] 


T Gi k Ti^{xLj , U/j) T Gj k T L (Ui, tL k ') Tl (Ui , 


d{uj) L x 

T L\ u ji u k) 

OX k 


d{uj)i 

dx k 


+ Cq £ 5 


ty 


(24) 


For VFDF2: 


d{uj) L _ 

d Xl 

d{uj) L d^iiUJ^h _ <HPh | _ ^LjuuUj) 

dt dxi dxj dx^Xi dx t 


(25) 
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d d $ 

— \T L (Ui,Uj)] + ■^■[{Uk) L T L (u l ,U j )\ = --^{T L (u u U } ,U k )} 

+ G ik T L (Uj,U k ) + G jk T L {Ui,U k ) - T L (Ui,U k ) ~ T L (Uj,U k ) + G 0 £ Sij. 


(26) 


It may be seen that the zeroth and first moment equations are identical (and exact); 
whereas the second central moment equations differ by the additional viscous term 
in VFDF1 (Eq. (24)). A comparison of these modeled equations with Eq. (5) shows 
that the GLM model implies: 


I (£ 


u 


£ &ij ) G 1 Ld 


t l {u u u 3 ) - -k 5 ij 


Ci — 1 + 2 ^°' 


(27) 


This is the same as the Rotta 33 model as shown by Pope. 34 There are two model 
constants in the VFDF equation. In RAS, typically 34,35 C e ~ 1, and Co ~ 2.\{C\ = 
4.15). As shown in Refs., 34,27 boundedness of the GLM coefficients C 0 > 0 guarantees 
that the SGS stress is realizable. 


4 EQUIVALENT STOCHASTIC SYSTEMS 

The solution of the VFDF transport equation provides all the statistical information 
pertaining to the velocity vector. The most convenient means of solving this equation 
is via the Lagrangian Monte Carlo scheme. The basis of this scheme relies upon 
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the principle of equivalent systems. 26,32 Two systems with different instantaneous 
behaviors may have identical statistics and satisfy the same PDF transport equation. 
In this context, the general diffusion process is considered via the following system of 
stochastic differential equations (SDEs), 26,36,37,31 


dXi{t) = Di(X(t),U(t)-,t) dt + B(X(t),U(t); t) dW t x {t), 
dUi{t) = dt + E{X(t)M{t)-,t) dW”(t) (28) 

+ F ij (X(t)Mt)\t) dW*(t), 


where X t and U x are probabilistic representations of the x and u , respectively. The 
coefficients D { and M, are the “drift” in the phase space of position and velocity, re- 
spectively. The terms B and E are the “diffusion” coefficients for physical and velocity 
spaces, respectively; and W* and W? denote independent Wiener-Levy processes. 38 
The tensor Fij represents the dependency between the velocity and physical spaces. 
This terms is needed to satisfy the Ito condition for B ^ 0. A comparison of the 
Fokker-Planck equation of Eq. (28) with the the modeled VFDF1 transport equation, 
Eq. (21) yields: 


Mi 


d(p) 


l , n . d 2 {ui) L 


dxi 


2 u 


dx k dx k 


+ Gij (Uj — (uj)i,)> Di =Ui, 


B = y/2^, E= F l} = \Flv 


i L 


di^i) 

dxj 


(29) 
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Therefore, the proper SDEs which represent VFDF1 in the Lagrangian sense are: 


dXi(t) = Ui(t) dt + dW*(t), 


dU.it) = 




dx. 


dx k dx k 


dt+ s/C^ dW?(t) 


+ sFu, dW ?(t). 


(30) 


This stochastic system is the same as that developed by Dreeben and Pope 29 31 for 
RAS. 

For VFDF2, due to the absence of diffusion in physical space we must have B = 0. 
Therefore, the corresponding SDEs are: 


dXi{t)=Ui(t) dt 

d(p) 


dUi{t) = 


l . d(a lk ) L 


dxi 


+ 


dx k 




dt+ y/Co£dW t v (t). 


(31) 


This system is the same as that suggested by Pope 26 and Haworth and Pope 27 for 
RAS. 

The primary difference between the two formulations VFDF1 and VFDF2 is due to 
molecular effects in the spatial diffusion of the VFDF. This is explicitly included 
in the VFDF1 formulation and is also present in the corresponding second moment 
equation. This difference is expected to be important in flows where viscous effects 
are important; e.g. flow near solid boundaries. 29-31 Both of these formulation are 
considered in our numerical simulations as discussed below. 
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5 NUMERICAL SOLUTION PROCEDURE 


Numerical solution of the modeled VFDF transport equation is obtained by a La- 
grangian Monte Carlo procedure. The basis of this procedure is the same as that in 
RAS 39 " 41 and in previous LES/FDF. 6,8 But there are some subtle differences which 
are explained here. In the Lagrangian description, the VFDF is represented by an 
ensemble of N statistically identical Monte Carlo particles. Each of these particles car- 
ries information pertaining to its velocity U^ n ^(f) and position X^ n ^(t), n — 1, 2, ... N. 
This information is updated via temporal integration of Eq. (30). The simplest means 
of performing this integration is via the Euler-Maruyamma approximation 42 

X?(t k+ i) = X?(t k ) + D?(t k )At + B n (t k ){At) l/2 Q(t k ), 

U?(t k+ 1) = U?(t k ) + M?(t k )At + E n (t k )(At)^ 2 C(tk) ( 32 ) 

+ F”(t k )(At) l ' 2 q(t k ), 

where D?(t k ) = A(X (n) (^), U<") (<*);«), B n (t k ) = B(X^(t k ),U^(t k )-,t), ... and 
Wk),Q(t k ) are independent standardized Gaussian random variables. This formu- 
lation preserves the Markovian character of the diffusion processes 43,44 and facilitates 
affordable computations. Higher order numerical schemes for solving Eq. (30) are 
available, 42 but one must be cautious in using them for LES. 6 Since the diffusion 
term in Eq. (28) strongly depends on the stochastic processes, the numerical scheme 
must be consistent with Ito-Gikhman 45,46 calculus. Equation (32) exhibits this prop- 
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erty. 


The statistics are evaluated by consideration of the ensemble of particles in a “finite 
volume” centered at a spatial location. This ensemble provides “one-time” statistics. 
This finite volume is characterized by a cubic box of length A E . This is necessary 
as, with probability one, no particle will coincide with the point as considered. 32 
Here, a cubic box of size A# is used to construct the ensemble mean, variances and 
covariances of the velocity vector. These values are used in the finite difference LES 
solver of Eq. (4) as described below. 

The SGS dissipation rate and the SGS mixing frequency as required in the solution 
of the VFDF are evaluated on the finite difference grid points and interpolated to the 
particle’s location. Ideally, for reliable Eulerian statistics and minimum numerical 
dispersion, it is desired to have the size of the sample domain infinitesimally small 
(i.e. A e -» 0) and the number of particles within this domain infinitely large. That 
is, 


P L {v,x,t) 


i 

Ne~> oo 
A e — ^0 


V NE (v;x.,t) = <5(«-U (n) ), 

NE neA E 


(33) 


where Vn e is the Eulerian PDF constructed from the particle ensemble, n € As 
denotes the particles contained in an ensemble box of length A# centered at x; and 
N e is the total number of particles within the box. With a finite number of particles, 
obviously a larger A# is needed. This compromise between the statistical accuracy 
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and dispersive accuracy implies that the optimum magnitude of A^ cannot, in general, 
be specified a priori. 26,11 This does not diminish the capability of the procedure, but 
exemplifies the importance of the parameters governing the statistics. 

To provide an estimate of the proper A# size, a “point estimator” procedure is con- 
sidered. With this procedure, the mean values (the first moments of the VFDF) are 
evaluated by ensemble averaging, and spatial variations of these mean values within 
the box are ignored. With the discrete representation (Eq. (32)), the first two mo- 
ments in this procedure are evaluated via: 


(w*)l * 


Nj? — > oo 
A e — ^0 


N e 


Y ^ (n) = 


ne Af 


T L (Ui,Uj) 


< 

Ne~+ °° 

Ag->0 


Y, W |n) - (Uih) 

e ne Ae 


(34) 


The point estimator is obviously subject to both statistical errors and dispersive errors 
for A e / 0. 

To determine the pressure field, the “mean field solver” is based on the “compact pa- 
rameter” finite difference scheme of Carpenter. 47 This is a variant of the McCormack 48 
scheme in which fourth-order compact differences are used to approximate the spatial 
derivatives, and a second order symmetric predictor-corrector sequence is employed 
for time discretization. The numerical algorithm is a hyperbolic solver which consid- 
ers a fully compressible flow. Here, the simulations are conducted with a low Mach 
number (M < 0.3) to minimize compressibility effects. All the finite difference opera- 
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tions are conducted on fixed and equally sized grid points. The transfer of information 
from these points to the location of the Lagrangian particles is conducted via interpo- 
lation. A second-order (bilinear) interpolation scheme is used for this purpose. The 
results of previous work indicate no significant improvements with the use of higher 
order interpolation schemes. 6 

The mean-field solver also determines the filtered velocity field. That is, there is a 
“redundancy” in the determination of the first filtered moments as both the finite 
difference and the Monte Carlo procedures provides the solution of this field. This 
redundancy is actually very useful in monitoring the accuracy of the simulated results. 
Detailed discussions pertaining to this issue are provided in Refs. 8,39 41 

To establish the consistency of the VFDF solver, another LES is also conducted 
in which the modeled transport equations for the filtered velocity and the general- 
ized SGS stresses are solved strictly via the finite difference scheme. These simula- 
tions are referred to as LES-FD and are only applied for the case corresponding to 
VFDF2. That is, Eqs. (25) and (26) are considered. Since the SGS transport terms 
r L (ui,Uj,Uk ) are unclosed in Eq. (26), the values corresponding to these terms are 
taken from the Monte Carlo solver and substituted in the SGS stress transport equa- 
tions. The attributes of all of the scheme are summarized in Table 1, with further 
discussions in Refs. 6,39-41 
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6 RESULTS 


A Flows Simulated 

Simulations are conducted of a two-dimensional (2D) planar jet, and a 3D tempo- 
rally developing mixing layer. The jet flow simulations are conducted primarily for 
establishing the consistency of the Lagrangian Monte Carlo solver. For this purpose, 
2D simulations are sufficient. To analyze the overall performance of the VFDF and 
to demonstrate its full capabilities and drawbacks, 3D simulations are required. 

In the planar jet, a fluid issues from a jet of width D into a co-flowing stream with 
a lower velocity. The size of the domain in the streamwise (x) and cross-stream 
( y ) directions are 0 < x < 14 D and —3.5 D < y < 3.5D. The ratio of the co- 
flowing stream velocity to that of the jet at the inlet is kept fixed at 0.5. A double- 
hyperbolic tangent profile is utilized to assign the velocity distribution at the inlet 
plane. The formation of the large scale coherent structures are expedited by imposing 
low amplitude perturbations at the inlet. In the finite difference simulations, the 
characteristic boundary condition procedure of Ref. 49 is used at the inlet, free-shear 
boundary conditions are used at the free-streams and the pressure boundary condition 
of Ref. 50 is used at the outflow. 

The temporal mixing layer consists of two parallel streams traveling in opposite direc- 
tions with the same speed. 51-53 A hyperbolic tangent profile is utilized to assign the 
velocity distribution at the initial time. The simulations are conducted for a cubic 
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box, 0 < x < L, -L/2 < y < L/2, 0 < z < L, where x, y and z denote the stream- 
wise, the cross-stream and the spanwise directions, respectively; and the length, L 
is specified such that L = 2 Np X u , where N P is the desired number of successive vor- 
tex pairings and A u is the wavelength of the most unstable mode corresponding to 
the mean streamwise velocity profile imposed at the initial time. The flowfield is 
parameterized in a procedure somewhat similar to that by Vreman et al. 20 The for- 
mation of the large scale structures are expedited through eigenfunction based initial 
perturbations. 54,55 This includes two-dimensional 52,20,56 and three-dimensional 52,57 
perturbations with a random phase shift between the 3 D modes. This results in the 
formation of two successive vortex pairings and strong three-dimensionality. 

The flow variables are normalized with respect to selected reference quantities. In 
the jet flow, the jet exit velocity and the jet width are the reference scales. In 
the temporal mixing layer, the reference length is the half initial vorticity thickness, 
L r = (A = , where ( u x ) L is the Reynolds averaged value of the 

2 ' \d{u\)l/dy\max 

filtered streamwise velocity and A U is the velocity difference across the layer). The 
reference velocity is U r = A U / 2 . 

B Numerical Specifications 

All finite difference simulations are conducted on equally-spaced grid points with grid 
spacings Ax = A y = Az(for 3D) = A. The resolution for LES of the planar jet 
consists of 201 x 101 grid points. This allows simulations with a Reynolds number 
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Re = = 14, 000. The simulations of the temporal mixing layer are conducted 

on 193 3 and 33 3 points for DNS and LES, respectively. This allows simulations with 
Re = 2d*. = 50. 

V 

To filter the DNS data, a top-hat function 21 of the form below is used 


N D 


G(x' - x) = G(x- - Xi) 


i=l 


G(Xi - x^ = < 


ST 


0 lx! - xA > 


(35) 


in which No denotes the number of dimensions, and A*, = 2A. 58 No attempt is made 
to investigate the sensitivity of the results to the filter function 24 or the size of the 
filter. 59 

For VFDF simulations of the temporal mixing layer, the Monte Carlo particles are 
initially distributed throughout the computational region. For the jet flow, the par- 
ticles are supplied in the inlet region — 1.75D < y < 1.75D. As the particles con- 
vect downstream, this zone distorts as it conforms to the flow as determined by the 
hydrodynamic field. The simulation results are monitored to ensure the particles 
fully encompass and extend well beyond regions of non-zero vorticity with an ap- 
proximately uniform particle number density. All simulations are performed with a 
uniform “weight” 26 of the Monte Carlo particles. In the temporal mixing layer, due 
to flow periodicity in the streamwise and spanwise directions, if the particle leaves 
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the domain at one of these boundaries new particles are introduced at the other 
boundary with the same compositional values. In the cross-stream directions, the 
free-slip boundary condition is satisfied by the mirror-reflection of the particles leav- 
ing through these boundaries. In the planar jet, new particles are introduced through 
the inlet boundary at a rate proportional to the local flow velocity and with a velocity 
makeup dependent on the cross-stream direction only. When the particles leave the 
computational domain at the outflow, they are no longer tracked. The density of the 
Monte Carlo particles is determined by the average number of particles Ne within 
the ensemble domain of size A e x A^(x A#). The effects of both of these parame- 
ters are assessed to ensure the consistency and the statistical accuracy of the VFDF 
simulations. 

All results are analyzed both “instantaneously” and “statistically.” In the former, 
the instantaneous contours (snap-shots) and scatter plots of the variables of interest 
are analyzed. In the latter, the “Reynolds-averaged” statistics constructed from the 
instantaneous data are considered. In the spatially developing flows this averaging 
procedure is conducted via sampling in time. In the temporal mixing layer, the 
statistics are constructed by spatial averaging over the x - z plane of statistical 
homogeneity. All Reynolds averaged results are denoted by an overbar. 
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C Consistency and Convergence Assessments 


The objective of this section is to demonstrate the consistency of the VFDF formula- 
tion and the convergence of its Monte Carlo simulation procedure. For this purpose, 
the results via VFDF and LES-FD are compared against each other. Since the accu- 
racy of the finite difference procedure is well-established (at least for the first order 
filtered quantities), such a comparative assessment provides a good means of assessing 
the performance of the Monte Carlo solution of the VFDF. This assessment is done 
for 2D simulations and for VFDF2 as they are sufficient for establishing consistency 
and convergence. To do so, the statistical results obtained from the Monte Carlo 
simulations of Eq. (31) are compared with the finite difference solution of Eqs. (25) 
& (26). Also, no attempt is made to determine the appropriate values of the model 
constants; the values suggested in the literature are adopted 34 Co = 2.1(Ci = 4.15) 
and C £ = 1. 

In Fig. 1, the instantaneous contour plots of the vorticity are shown as determined by 
(a) VFDF2 and (b) LES-FD. This figure provides a simple visual demonstration of 
the consistency of the VFDF2. Scatter plots of ( u)l vs. ( v)i are presented in Fig. 2. 
The correlation and regression coefficients (denoted, respectively, by p and r on these 
figures) are insensitive to A e- Figures 3 and 4 show the Reynolds averaged values of 
the streamwise velocity and several components of the SGS stress tensor for several 
values of A E , with Ne = 40 kept fixed. It is observed that the first filtered moments 
as obtained by VFDF agree very well with those via LES-FD even for large A e values. 
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However, smaller A# values are required for convergence of the VFDF predicted SGS 
stresses to those by LES-FD. The relative difference between the L 2 norms of all of 
the components of the SGS tensor as a function of (A#/A) 2 is presented in Fig. 5. 
Extrapolation to A^ = 0 shows that the “error” goes to zero as A e 0. 

The influence of Ne on the first two moments is shown in Figs. 6 and 7. It is observed 
that Ne does not have a significant influence on the first moments, but does slightly 
influence the second moments. In all the cases considered, Ne > 40 yields reliable 
predictions, consistent with previous consistency and convergence assessments of the 
scalar FDF. 6,8 

All the subsequent simulations are conducted with A e = A/2 and Ne = 40. 

D Comparative Assessments of the VFDF 

The objective of this section is to analyze some of the characteristics of the VFDF via 
comparative assessments against DNS data. This assessment is done via both a priori 
and a posteriori analyses. In the former, the DNS results are used to determine the 
range of the empirical constants appearing in the VFDF sub-closures. In the latter, 
the final results as predicted by the VFDF are directly compared with those obtained 
by DNS. The procedure is similar to that in Ref. 20 and considers the 3D temporal 
mixing layer. 

In addition to VFDF, three other LES are conducted with (1) no SGS model, (2) the 
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Smagorinsky 16 ’ 60 SGS closure, and (3) the dynamic Smagorinsky 17-19 model. In the 
case with no model, the contribution of the SGS is completely ignored, i.e. Uj ) = 

0. In this case, the numerical errors amount to an implied model. But as indicated 
in Ref., 20 this case is included to provide a point of reference for the other closures. 
The Smagorinsky model is, 16 ’ 61 


tl{u u Uj) - -k 8ij = -2 v t Sij , 

1 / d(ui) L d(uj) L \ 
lj ~ 2 V dxj dxi ) 

Vt = a a 2 5. 


(36) 


C v = y/2 0.17 2 « 0.04, S = s/S lJ S l] and A L is the characteristic length of the 
filter. This model considers the anisotropic part of the SGS stress tensor Oy = 
TL(u{,Uj ) — 2/3 k 5 tJ . The isotropic components are absorbed in the pressure field. 
The dynamic version of the Smagorinsky model provides a means of approximating 
C v as suggested in Refs. 17 19 The procedure for the implementation of this model in 
the 3D temporal mixing layer LES is described by Vreman; 20 thus it is not repeated 
here. (See Refs. 11,23,62,63 for recent reviews on SGS closure strategy.) 

In addition to the resolved velocity field, the primary integral statistical quantities 
considered for comparative assessments are: 


Ef 



(37) 
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( 38 ) 


Pk = 

J 

1 p k dx, with p k 

E v = 

j e u dx, with e u 

B k = 

[ min(0,p k ) dx. 


d{ui) 


i/L 


d(uih ( 9(uj )_l 


+ 


d{u k ) L \ 
dxi ) ’ 


Ef is the kinetic energy of the resolved field, t v represents the viscous molecular dis- 
sipation rate directly from the filtered field, P k is the production rate of the SGS 
kinetic energy (or the rate of energy transfer from the resolved filtered motion to 
the SGS motion), and B k is the total backscatter. 64-66 The resolved molecular dis- 
sipation rate is always positive (by definition), but the production rate of the SGS 
kinetic energy can be locally negative. This backscatter is not represented in the 
Smagorinsky model. The dynamic model is potentially capable of accounting for 
it, but at the expense of causing numerical instabilities. In the implementation 
of the dynamic model used here, backscatter is avoided by averaging the numera- 
tor and denominator of the expression determining C„ (Refs. 19,20 ) over the homoge- 
neous directions. If negative values are still present, they are set equal to zero. 20,63 
The “resolved” components of the Reynolds-averaged stress tensor are denoted by 
R~j where R tj = ((Ui) L - Wh)(( U j)L ~ Wjh)- The “total” Reynolds stresses 
are denoted by f~ where r xj = (U t - Ul){Uj - TJ]). These are approximated by 
jrr ~ r~. + riJuuUj ). 20, 67 ' 68 In DNS, the total stresses are evaluated directly and the 
results indicate that R~j + ^(uj.Uj) does indeed approximate r~ with a maximum 
error of less than 10% 
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Figure 8 shows the distribution of the particle number density within the whole 
computational domain. Assuring an approximately uniform distribution, the values 
of the moments within local ensembles are compared with those of filtered DNS data. 
These DNS data are transposed from the original high resolution 193 3 points to the 
low resolution of 33 3 points, and then are compared with LES results on these coarse 
points. 

The DNS data are also used to make a priori estimates of the model constants. The 
primary terms which require closure are the SGS dissipation and the velocity-pressure 
scrambling tensors. The model equation (Eq. (20)) involving C £ is in a scalar form. 
For an estimate of C x (thus C 0 ), we consider the following norm of the corresponding 
closure (Eq. (27)): 


|| - II ij - (Sij ~^£ <5q)|| ~ Ci u || r L {u t , Uj) - | k <y|, (39) 

where ||WA,|| = s/W^Wij- To estimate the coefficients, a linear regression is performed 
on all the data points at each computational time step. The optimized constants as 
obtained in this way are denoted by C e and C x . This procedure is also followed for 
the Reynolds averaged data, with the optimized models obtained in this way denoted 
by C e and C x . The temporal variations of these estimated values are shown in Fig. 9. 
The non-uniformity of the coefficients indicates the “non-universality” of the models. 
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This is expected as the flow evolves from an initially smooth laminar state to a strong 
three-dimensional state (at t « 40) before the action of the small scales becomes 
significant. The closures as adopted are not fully suitable for application in all of 
these flow regions. Nevertheless, Fig. 9 indicates that the values for these coefficients 
as suggested in RAS, i.e. C x « 4.15, C e « 1 are reasonable, at least within the 
turbulent regime. The influences of these parameters are further investigated via a 
posteriori analysis of the results as discussed below. 

Figures 10 and 11 show the contours of the spanwise and the streamwise components 
of the vorticity field, respectively, at time t = 80. By this time, the flow has gone 
through several pairings and exhibits strong 3D effects. This is evident by the for- 
mation of large scale spanwise rollers with presence of counter-rotating streamwise 
vortex pairs in all the simulations. The results via the no-model indicate too many 
small-scale structures which clearly are not captured accurately on the coarse grid. 
The amount of SGS diffusion with the Smagorinsky model is very significant at initial 
times. Due to this dissipative characteristics of the model, the predicted results are 
too smooth and only contain the large scale structures. The vortical structures as 
depicted by the dynamic Smagorinsky and the VFDF are very similar and predict 
the DNS results better than the other two models. The results obtained by VFDF1 
and VFDF2 are virtually indistinguishable from each other. This is expected, due to 
the lack of importance of molecular effects in this free shear flow. 

The Reynolds averaged values of the streamwise velocity and the temporal variations 
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of the momentum thickness: 


1 r L/ 2 

s m (t) = -l ( 1 - {u) L ) ( 1 + (u) L ) dy, (40) 

4 J-L/2 

are shown in in Figs. 12 and 13, respectively. In Fig. 12 the Reynolds averaged values 
of both filtered and unfiltered DNS data are considered and are shown to be essentially 
equivalent. Therefore, the latter are not shown in subsequent figures. The dissipative 
nature of the Smagorinsky model at initial times resulting in a slow growth of the 
layer is shown. Several values of the model parameters (Co, C e ) are considered in the 
VFDF simulations. It is observed that as the magnitude of C £ decreases, the initial 
rate of the layer’s spread is higher. With the exception of the case with C £ = 0.5 and 
the Smagorinsky model, all the other VFDF cases, the dynamic Smagorinsky and the 
no-model yield a similar rate of layer’s growth at late times. 

The temporal variations of the resolved kinetic energy and all of the terms defined in 
Eq. (38) are shown in Fig. 14. The overall features displayed in this figure are similar 
to those reported by Vreman et al. 20 for the no model, the Smagorinsky model and the 
dynamic Smagorinsky model. The initial rate of decay of the resolved kinetic energy 
for the Smagorinsky model is the highest. This is due to the excessive production 
of the SGS kinetic energy by this model in the transitional region, and explains 
the reason for the lack of small scales in the vortical structures as discussed before. 
For all the other models the initial rate of decrease of the resolved kinetic energy 
is small and increases as the flow develops. The trend portrayed by DNS results is 
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best captured by the VFDF simulations. For the no model case the only means of 
dissipation of the resolved kinetic energy is through molecular action and numerical 
dissipation which become significant at later stages due to presence of a large amount 
of small scales. In this case, the amount of numerical dissipation is the highest. For 
all the other closures, the productions rate of the SGS kinetic energy is larger than 
the molecular dissipation as the flow develops. The dynamic Smagorinsky and the 
no-model simulations predict the same initial rate of decay for the resolved kinetic 
energy. This is due to low initial values of P k predicted by the dynamic Smagorinsky 
model. After t = 40 the amount of P k as predicted by the dynamic model is more 
than that of molecular dissipation by the no-model. Thus the rate of decay of the 
resolved kinetic energy becomes higher for the dynamic model and is closer to that 
obtained by DNS. 

With the exception of the no-model case, all the simulations predict similar trends for 
the molecular dissipation. The magnitude of this dissipation as predicted by VFDF 
changes slightly with the variation of the model parameter. The production rate of the 
SGS kinetic energy depends more strongly on the model coefficients; as C e decreases, 
the peak magnitude of P k is larger. The Smagorinsky model does not adequately 
predict P k , and the dynamic model yields better predictions at long times. The overall 
trends are best predicted by VFDF. The same is true in capturing the backscatter 
phenomenon. By design, the backscatter is identically zero in the Smagorinsky and 
the dynamic Smagorinsky model. But VFDF is capable of capturing it, and its extent 
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is controlled by the model parameters. In this regard it is important to note that 
there are no numerical instability problems in the VFDF solver for negative £* values. 
However, the amount of predicted backscatter is less than that of DNS and its relative 
magnitude is less than those of P* and E v . 

Several components of the planar averaged values of the SGS anisotropy tensor, = 
T h (u u Uj) - 2/3 k 5ij are presented in Figs. 15 & 16. Both the Smagorinsky and 
the dynamic model under-predict all of the components of this stress. The VFDF 
predictions are much more satisfactory. In this regard, the VFDF is expected to be 
more effective than the other closures for LES of reacting flows since the extent of 
SGS mixing is heavily influenced by SGS convection. 69 ’ 70 “Optimum” values for C e 
and Co cannot be suggested to predict all of the components of this tensor at all 
times, but it is obvious that there is too much SGS energy with C e — 0.5. 

Several components of the resolved stress tensor Rij are shown in Figs. 17 & 18. As 
expected, the performance of the Smagorinsky model is not very good as it does not 
predict the spread and the peak value of the resolved Reynolds stresses. None of 
the other models show a distinct superiority in predicting the DNS results. The no- 
model and the dynamic Smagorinsky model predict large peak values at the middle 
of the layer. The VFDF predicts both the spread and the peak values reasonably 
well. The results for small C e values are not shown since the amount of energy in 
the resolved scale decreases too much in favor of the increase of the SGS stress (as 
shown in Figs. 15 & 16). The cross-stream variations of the total Reynolds stress r [2 
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are presented in Figs. 19. The peak values by the no-model simulations are again 
the highest. The dynamic model and VFDF perform similarly and capture the DNS 
trends equally well. 

E Comparison with Previous Investigations 

All of the results obtained here by DNS, and LES via the Smagorinsky and the dy- 
namic Smagorinsky models agree very well with those of Vreman et al. 20 The slight 
differences are due to the non-identical flow initializations, and the different compu- 
tational methodologies employed in the two simulations. To compare with results 
of other investigations, simulations are conducted of another temporally developing 
mixing layer with Re = 500 in a larger computational domain, L r = 120. An ini- 
tial forcing of the form Ae~^ v ^ 2 is used, where A is a uniformly distributed random 
number with an amplitude of 0.05. Rogers and Moser 60 perform DNS of a high Re 
number flow on 512 x 210 x 192 spectral points. The results of these simulations 
are in excellent agreements with laboratory data of Bell and Mehta. 71 Here, LES is 
conducted of this flow via the dynamic Smagorinsky model. 

The profiles of the mean streamwise velocity and several components of the resolved 
stresses at t = 250 are presented in Figs. 20 and 21, respectively. In these figures, 
£ = y/S m (t) and the symbols denote the experimental data 71 at several streamwise 
locations. The good agreement with these data also indicates good agreement with 
DNS results of Rogers and Moser. 72 
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F Computational Requirements 


The total computational times associated with simulations of the 3D temporal mixing 
layer are shown in Table 2. Expectedly, the overhead associated with the VFDF 
simulation is extensive as compared to the other models; nevertheless this requirement 
is significantly less that of DNS. This overhead was tolerated in present simulations, 
but can be reduced with utilization of an optimum parallel simulation procedure. 
This has been discussed for use in PDF 73 and is recommended for future VFDF 
simulations. 


7 SUMMARY AND CONCLUDING REMARKS 

The filtered density function (FDF) methodology 1 has proven very effective for large 
eddy simulation (LES) of turbulent reacting flows. 3,6-11 In all of previous contribu- 
tions, the LES/FDF of only the scalar quantities are considered. The objective of the 
present work is to develop the FDF methodology for LES of the velocity field. For 
this purpose, a methodology termed the velocity filtered density function (VFDF) 
is developed. The VFDF is basically the probability function (PDF) of the subgrid 
scale (SGS) velocity vector. The exact transport equation governing the evolution of 
the VFDF is derived. It is shown that the effects of SGS convection in this equation 
appears in a closed form. The unclosed terms in this transport equation are modeled 
via two formulations: VFDF1 and VFDF2. The primary difference between the two 
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models is the inclusion of the molecular diffusion in the spatial transport of the VFDF 
in the first formulation. The closure strategy in the formulation similar to that in 
PDF methods in Reynolds averaged simulation (RAS) procedures . 32 In this way, the 
VFDF formulation is at least equivalent to a second-order moment SGS closure. 

The modeled VFDF transport equations are solved numerically via a Lagrangian 
Monte Carlo scheme in which the solutions of the equivalent stochastic differential 
equations (SDEs) are obtained. Two Monte Carlo procedures are considered. The 
schemes preserve the Ito-Gikhman nature of the SDEs and provide a reliable solution 
for the VFDF. The consistency of the VFDF formulation and the convergence of 
its Monte Carlo solutions are assessed. This is done via comparisons between the 
results obtained by the Monte Carlo procedure and the finite difference solution of 
the transport equations of the first two filtered moments of VFDF (LES-FD). With 
inclusion of the third moments from the VFDF into the LES-FD, the consistency and 
convergence of the Monte Carlo solution is demonstrated by good agreements of the 
first two SGS moments with those obtained by LES-FD. 

The VFDF predictions are compared with those with LES results with no SGS model, 
with the Smagorinsky 16 SGS closure, and with the dynamic Smagorinsky 17-19 model. 
All of these results are also compared with direct numerical simulation (DNS) results 
of a three-dimensional, temporally developing mixing layer in a context similar to 
that conducted by Vreman et al . 20 This comparison provides a means of examining 
some of the trends and overall characteristics as predicted by LES. It is shown that 
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the VFDF performs well in predicting some of the phenomena pertaining to the SGS 
transport. The magnitude of the SGS Reynolds stresses as predicted by VFDF is 
significantly larger than those predicted by the other SGS models and much closer 
to the filtered DNS results. The temporal evolution of the production rate of the 
SGS kinetic energy is predicted well by VFDF as compared with those via the other 
closures. The VFDF is also capable of accounting the SGS backscatter without any 
numerical instability problems, although the level predicted is substantially less than 
that observed in DNS. 

The results of a priori assessment against DNS data indicates that the values of the 
model coefficients as employed in VFDF (Co and C £ ) are of the range suggested in the 
equivalent models previously used in RAS. The results of a posteriori assessments via 
comparison with DNS data does not give any compelling reasons to use values other 
than those suggested in RAS, C 0 = 2.1, C e = 1. However, small values of C E are not 
acceptable as they would yield too much of SGS energy relative to that within the 
resolved scales. 

Most of the overall flow features, including the mean velocity field and the resolved 
and total Reynolds stresses as predicted by VFDF are similar to those obtained via the 
dynamic Smagorinsky model. This is interesting in view of the fact that the model 
coefficients in VFDF are kept fixed. It may be possible to improve the predictive 
capabilities of the VFDF by two ways: (1) development of a dynamic procedure to 
determine the model coefficients, and/or (2) implementation of higher order closures 
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for the generalized Langevin model parameter Gi } (see Ref. 34 ). 

Work is in progress towards developments of a joint velocity-scalar FDF for LES 
of reacting flows. Compared to standard LES, this approach has the advantage of 
treating reaction in a closed form; and, compared to scalar FDF 6,8 has the advantage 
of treating convective transport (of momentum and species) in closed form. These 
modeling advantages have an associated computational penalty. For the cases consid- 
ered here, VFDF is more expensive computationally than the dynamic Smagorinsky 
model by a factor of 15. It is expected that VFDF will not be more expensive than 
scalar FDF, at least for reacting flows with many species. 
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Table 1: Recapitulation of the VFDF solution procedures. 


49 












Resolution 

n e 

Normalized CPU time 

DNS 

193 

X 

193 

X 

193 

- 

178 

VFDF1 

33 

X 

33 

X 

33 

40 

33.6 

VFDF2 

33 

X 

33 

X 

33 

40 

30 

Dynamic Smagorinsky 

33 

X 

33 

X 

33 

- 

2.19 

Smagorinsky 

33 

X 

33 

X 

33 

- 

1.05 

No model 

33 

X 

33 

X 

33 

- 

1 


Table 2: Computer requirements for the 3D temporal mixing layer. 1 unit corresponds 
to 1657.2 seconds of CPU time on the SGI Origin 2000. 
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FIGURE CAPTIONS 


Figure 1. Plot of the vorticity field contours, (a) VFDF2, (b) LES-FD. A e — 
A, N e = 40. 

Figure 2. Scatter plots of the filtered velocity field as obtained via VFDF2 vs. 
LES-FD. (a), (u) L ; (b), (v) L . A E = A , N E = 40. 

Figure 3. Reynolds averaged values of the filtered streamwise velocity, (a) 
Cross-stream variations at x = 7, (b) streamwise variation at y = 0 (center- 
line). N e = 40. 

Figure 4. Cross-stream variations of the Reynolds averaged values of some of 
the components of the SGS stress tensor at x = 7 with N e = 40. The LES-FD 
results are obtained with A# = 0.5A, N e = 40. 

Figure 5. Percentage of the relative difference between the L 2 norms of the 
stresses as a function of (a) x = 2.8, (b) x — 7, (c) x — 11.2. 

Figure 6. Cross-stream variations of the Reynolds averaged values of the filtered 
streamwise velocity at x = 7. The LES-FD results are obtained with A# = 
0.5A, N e = 40. 

Figure 7. Cross-stream variations of the Reynolds averaged values of some of 
the components of the SGS stress tensor at x = 7. The LES-FD results are 
obtained with A e = 0.5A, N e = 40. 
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Figure 8. Particle number density in VFDF2 simulation at t = 60. The iso- 
surface corresponds to Ng — 40 set as initial conditions. Co = 2.1, C £ = 1. 

Figure 9. Time variation of the model coefficients as obtained from a priori 
analysis of the DNS data. 

Figure 10. Contour plots of the spanwise component of the vorticity at z = 
0.75 L/L r , t = 80. (a) filtered DNS, (b) no model, (c), Smagorinsky model, 
(d) dynamic Smagorinsky model, (e) VFDF2, Co = 2.1, C £ = 1, (f) VFDF1, 
C 0 = 2.1, C e = 1. 

Figure 11. Contour plots of the streamwise component of the vorticity vector at 
x = 0.25 L/L r , t = 80. (a) filtered DNS, (b) no model, (c) Smagorinsky model, 
(d) dynamic Smagorinsky model, (e) VFDF2, Co = 2.1, C £ = 1, (f) VFDF1, 
Co = 2.1,C e = l. 

Figure 12. Cross-stream variations of the Reynolds averaged values of the 
streamwise velocity at t = 70. 

Figure 13. Temporal variations of the momentum thickness. 

Figure 14. Temporal variations of (a) total resolved kinetic energy, (b) SGS 
kinetic energy production rate, (c) total backscatter, (d) total resolved dissipa- 
tion. 

Figure 15. Cross stream variations of some of the components of Ojj at t = 60. 
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Figure 16. Cross-stream variations of some of the components of ai] at t = 80. 

Figure 17. Cross-stream variations of some of the components of R tJ at t = 60. 

Figure 18. Cross-stream variations of some of the components of R t j at t = 80. 

Figure 19. Cross stream variations of rj-J, (a) t = 60, (b) t = 80. 

Figure 20. Cross stream variation of the Reynolds averaged values of the stream- 
wise velocity at t = 250. Solid line denotes model predictions via the dynamic 
Smagorinsky model. Symbols denote experimental data of Bell and Mehta. 71 

Figure 21. Cross stream variations of the Reynolds averaged values of the 
streamwise velocity at t = 250. Solid lines denote model predictions via the 
dynamic Smagorinsky model. Symbols denote experimental data of Bell and 
Mehta. 71 
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(b) LES-FD 
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LES-FD 

Figure 2: Scatter plots of the filtered velocity field as obtained via VFDF2 vs. LES 
FD. (a), (u) L \ (b), (v) L . A e = A, N E = 40. (Gicquel et al., Physics of Fluids.) 
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Figure 3: Reynolds averaged values of the filtered streamwise velocity, (a) Cross- 
stream variations at x = 7, (b) streamwise variation at y — 0 (center-line). N E = 40 
(Gicquel et al., Physics of Fluids.) 
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Figure 4: Cross-stream variations of the Reynolds averaged values of some of the 
components of the SGS stress tensor at x = 7 with N E — 40. The LES-FD results 
are obtained with A E = 0.5A, N E = 40. (Gicquel et ah, Physics of Fluids.) 
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Figure 6: Cross-stream variations of the Reynolds averaged values of the filtered 
streamwise velocity at x = 7. The LES-FD results are obtained with A# = 
0.5A, N e = 40. 
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Figure 7: Cross-stream variations of the Reynolds averaged values of some of the 
components of the SGS stress tensor at x = 7. The LES-FD results are obtained 
with A e = 0.5A, N e = 40. (Gicquel et al., Physics of Fluids.) 


60 






Figure 8: Particle number density in VFDF2 simulation at t = 60. The iso-surface 
corresponds to Ne — 40 set as initial conditions. Co = 2.1, C e = 1. (Gicquel et al., 
Physics of Fluids.) 


61 



e variation of the model co 
ta. (Gicquel et al., Physics 






Figure 10: Contour plots of the spanwise component of the vorticity at z = 

0.75 L/L r , t = 80. (a) filtered DNS, (b) no model, (c), Smagorinsky model, (d) 

dynamic Smagorinsky model, (e) VFDF2, Co = 2.1, C e = 1, (f) VFDF1, Co — 2.1, 
C £ = 1. (Gicquel et ah, Physics of Fluids.) 
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Figure 11: Contour plots of the streamwise component of the vorticity vector at 
x = 0.25 L/L r , t = 80. (a) filtered DNS, (b) no model, (c) Smagorinsky model, (d) 
dynamic Smagorinsky model, (e) VFDF2, Co = 2.1, C £ — 1, (f) VFDF1, Co = 2.1, 
C e = 1. (Gicquel et al., Physics of Fluids.) 
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Figure 12: Cross-stream variations of the Reynolds averaged values of the streamwise 
velocity at t = 70. (Gicquel et ah, Physics of Fluids.) 
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Figure 13: Temporal variations of the momentum thickness. (Gicquel et al., Physics 
of Fluids.) 
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Figure 14: Temporal variations of (a) total resolved kinetic energy, (b) production 
rate of the SGS kinetic energy, (c) total backscatter, (d) total resolved dissipation. 
(Gicquel et al., Physics of Fluids.) 
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Figure 15: Cross stream variations of some of the components of a tJ at t — 60. 
(Gicquel et al., Physics of Fluids.) 
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Figure 18: Cross-stream variations of some of the components of R t j at t 
(Gicquel et al., Physics of Fluids.) 
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Figure 19: Cross-stream variations of r\ 2 - ( a ) t = 60, (b) t = 80. (Gicquel et al 
Physics of Fluids.) 
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Figure 20: Cross stream variation of the Reynolds averaged values of the streamwise 
velocity at t — 250. Solid line denotes model predictions via the dynamic Smagorin- 
sky model. Symbols denote experimental data of Bell and Mehta. 71 (Gicquel et al., 
Physics of Fluids.) 
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(a) 



(b) 




(d) 



Figure 21: Cross stream variation of the resolved Reynolds stresses at t — 250. Solid 
lines denote model predictions via the dynamic Smagorinsky model. Symbols denote 
experimental data of Bell and Mehta. 71 (Gicquel et al., Physics of Fluids.) 
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Abstract 

A methodology termed the “velocity-scalar filtered density function” (VSFDF) is developed and 
implemented for large eddy simulation (LES) of turbulent flows. In this methodology, the effects 
of the unresolved subgrid scales (SGS) are taken into account by considering the joint probability 
density function (PDF) of the velocity-scalar field. An exact transport equation is derived for the 
VSFDF in which the effects of the SGS convection and chemical reaction are closed. The unclosed 
terms in this equation are modeled. A system of stochastic differential equations (SDKs) which 
yields statistically equivalent results to the modeled VSFDF transport equation is constructed. 
These SDEs are solved numerically by a Lagrangian Monte Carlo procedure in which the Ito- 
Gikhman character of the SDEs is preserved. The consistency of the proposed SDEs and the 
convergence of the Monte Carlo solution are assessed by comparison with results obtained by a 
finite difference LES procedure in which the corresponding transport equations for the first two 
SGS moments are solved. The VSFDF results are compared with those obtained via other SGS 
closures, and all the results are assessed via comparison with data obtained by direct numerical 
simulation (DNS) of a temporally developing mixing layer involving transport of a passive scalar. It 
is shown that the values of both the SGS and the resolved components of all second order moments 
including the scalar fluxes are predicted well by VSFDF. The sensitivity of the model’s (empirical) 
constants are assessed and it is shown that the magnitudes of these constants are in t he same range 
as that typically employed in PDF methods. 

Corresponding Author: Peyman Givi, Tel: (716) 645-2593 (x. 2320). Fax: 716-645-3875. E-mail: 
givi eng.buffalo.edu. 

Keywords: large eddy simulation, filtered density function, probability density function, turbulent combus- 
tion, Monte Carlo simulation. 
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I. INTRODUCTION 


The probability density function (PDF) approach has proven useful for large eddy sim- 
ulation (LES) of turbulent reacting flows [1]. The formal means of conducting such LES is 
by consideration of the “filtered density function” (FDF) [2] which is essentially the filtered 
fine-grained PDF of the transport quantities. In all previous contributions, the “marginal” 
FDF of the scalars [3-15], or the marginal FDF of the velocity vector [16] are considered; 
see Givi [17] for a recent review. 

The objective of the present work is to extend the FDF methodology to account for the 
“joint” subgrid scale (SGS) velocity-scalar field. This is facilitated by consideration of the 
joint “velocity-scalar filtered density function” (VSFDF). With the definition of the VSFDF, 
the mathematical framework for its implementation in LES is established. A transport 
equation is developed for the VSFDF in which the effects of SGS convection and SGS 
chemical reaction (in a reacting flow) are closed. The unclosed terms in this equation are 
modeled in a fashion similar to those in the Reynolds-averaged simulation (RAS) procedures. 
A Lagrangian Monte Carlo procedure is developed and implemented for numerical simulation 
of the modeled VSFDF transport equation. The consistency of this procedure is assessed 
by comparing the first two moments of the VSFDF with those obtained by the Eulerian 
finite difference solutions of the same moments’ transport equations. The results of the 
VSFDF simulations are compared with those predicted by the Smagorinsky [18] closure. All 
the results are assessed via comparisons with direct numerical simulation (DNS) data of a 
three-dimensional (3D) temporally developing mixing layer involving transport of a passive 
scalar variable. 

II. FORMULATION 

For the general formulation, we consider an incompressible (unit density), isothermal, 
turbulent reacting flow involving N a species. The primary transport variables describing 
such a flow are the three components of the velocity vector iq(x,f) (i = 1,2,3), the pressure 
p(x, t), and the species’ mass fractions <^ Q (x, t) (a = 1, 2, . . . , A^). The equations which 
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govern the transport of these variables in space (ar,-) and time ( t ) are 


dui 

dt 
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dt 
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= 0 
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_ dp dcr tk 
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dxi dx k 
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where S a = S a ((f){x,t)) denotes the chemical reaction term for species o, and <j> = 
[<f> u <^> 2 -. ■ • • , denotes the scalar variable array. For an incompressible, Newtonian fluid, 
with Fick’s law of diffusion, the viscous stress tensor a lk and the scalar flux Jg are repre- 
sented by 


®ik 


V 




■n = 


r d<f> a 

dx k 


(2a) 

(2b) 


where v is the fluid kinematic viscosity and T = is the diffusion coefficient of all species 
with Sc denoting the molecular Schmidt number. We assume a constant value for u = T; 
that is Sc = 1. In reactive flows, molecular processes are much more complicated than 
portrayed by Eq. (2). Since the molecular diffusion is typically less important than that of 
SGS, this simple model is adopted with justifications and caveats given by in Refs. [19-21], 
Large-eddy simulation involves the spatial filtering operation [22-25] 

/ +oo 

f{x',t)Q{x\x)dx' (3) 

-OG 

where Q(x',x) denotes a filter function, and ( f(x,t )) is the filtered value of the transport 
variable f(x,t). We will consider a filter function that is spatially and temporally invariant 
and localized, thus: Q(x',x) = G(x' — x) with the properties G(x) > 0, G(x)dxl, and 
moments x m G(x)dx exist for m > 0. Applying the filtering operation to Eqs. (1) yields 


d (uk) 
dx k 

<9(u,) d(u k ) (uj) 
dt dxk 

d_(4>a) d (Uk) (4> a ) 
dt dx , t 


= 0 

d{p) y d 2 jm) _ dT(u k ,Uj) 
dx i dxkdxk dx/, 

d 2 (<j> Q ) dT(u k ,<f> a ) 
v dx k dx k dx k 


(4a) 

(4b) 

(4c) 


4 



where the second-order SGS correlations 


r(a,b) = ( ab ) - (a) ( b ) 


( 5 ) 


are governed by 
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In this equation, the third order correlations 

r(a,b,c) = (abc) — (a) r(b,c) 

- ( b ) r(a, c) - ( c ) r(a , 6) - (a) (6) (c) (7) 

are unclosed along with the other terms within square brackets. 


III. VELOCITY-SCALAR FILTERED DENSITY FUNCTION (VSFDF) 
A. Definitions 


The “velocity-scalar filtered density function” (VSFDF), denoted by P , is formally defined 
as 

/ +oo 

g{v,ijr, u(x', t), t)) G(x' - x)dx' (8) 

■oo 

3 N s 
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a = l 
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where 8 denotes the delta function, and are the velocity vector and the scalar array 
in the sample space. The term g is the “fine-grained” density [20, 26], hence Eq. (8) 
defines VSFDF as the spatially filtered value of the fine-grained density. With the condition 
of a positive filter kernel [27], P has all of the properties of the PDF [20]. For further 
developments it is useful to define the “conditional filtered value” of the variable Q(x,t) as 



u(x y t) = l 


v>) = (q|v,v>) 

P(v,ip;xJ) 


x) dx' 


( 10 ) 


Equation (10) implies the followings: 


(i) for Q(x, t) = c, 

(ii) for Q(x, t) = Q{u{x, t), </>(x, t)) 

(iii) Integral properties: 


(<2(x,0 1 = c 

| v, ip^j = Q(v,ip) 


(<5(x,<)) 



(^Q{x,t) v,ip^ P( 


(11a) 

(lib) 


v,ip;x, t)dvdip 


(11c) 


From Eqs. (11) it follows that the filtered value of any function of the velocity and/or scalar 
variables is obtained by its integration over the velocity and scalar sample spaces 


(Q{x 


/ + OG 

' 

■oo 


Q(v , ip)P(v, ip; x, t)dvdip 


( 12 ) 


B. VSFDF Transport Equations 

To develop the VSFDF transport equation, we consider the time derivative of the fine- 
grained density function (Eq. (9)) 


dg _ _ / du^dp^ d<f> n dp \ 

dt \ dt dvk dt dipa) 

Substituting Eqs. (lb-lc), and Eqs. (2a-2b) into Eq. (13) we obtain 

dg dukp _ ( dp_ _ (Pug \ 9 q_ _ f d 2 <t>a s , , A dp 
dt dxk \dxi V dxkdxk) dv { \ dx k dxk ° / dip 0 


(13) 


(14) 
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Integration of this according to Eq. (8), while employing Eq. (10) results in 


c)P_ dvkP_ 
dt + dxk 


+ 


d (g) dP d 
dxk dv k di’ a 

d_ 
dv k 
d_ 

dvi 

d 


[SM)P] 


dtl> a 


dp 

dx k 

d 2 Ui 

dxkdxk 
v d^a 
dxkdxk 


v,ip) - 


dx k 
\ 
P 


v,1p) P 


( 15 ) 


This is an exact transport equation for the VSFDF. It is observed that the effects of convec- 
tion (second term on LHS) and chemical reaction (the last term on RHS) appear in closed 
forms. The unclosed terms denote convective effects in the velocity-scalar sample space. 
Alternatively, the VSFDF equation can be expressed as 


dP_ dvkP 
dt dxk 


d 2 P d ( p ) dP d 
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dxkdxk dxk dvk 
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dvk 

d 2 


dp 


v,ip 


d^ a 

d( P ) 


[SMP] 


dvidvj 
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d 2 


dx k 

dui duj 
dx k dx k 
dui dip, 


dx k 

v, ijj ) P 


dvidtfra [ 

d 2 


dxk dxk 
dipc drpp 


d^ad^p \\ dx k dx k 


v,1p) P 
P 


(16) 


This is also an exact equation, but the unclosed terms are exhibited by the conditional 
filtered values of the dissipation fields as shown by the last three terms on the RHS. 
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C. Modeled VSFDF Transport Equation 


For closure for the VSFDF transport equation, we consider the general diffusion process 
[28], given by the system of stochastic differential equations (SDEs): 

dX?(t) = Df(X + ,\J + ,cj ) + -t)dt + Bf J (X + ,V + ,<f> + ] t)dWf(t) 

+ F* U (X + , U + , t)dW?{t) + Ft* (X + , U + , </> + ; t)dWf(t) (17a) 
dU?{t) = Df(X+,U + ,^ + ;0^ + ^(X + ,U + ,^ + ;t)dHf(<) 

+ F^ x (X\l] + ,<l ) + G)dWf(t) + F^(X\lJ + ^ + G)dWf(l) (17b) 

d<j>l{t) = Di(X + ,\J + ,<f>+;t)dt + B+ j (X+,U + ,4> + -,t)dWj‘(t) 

+ (X + , U + , <f> + ] t)dW* (t) + F^{X + ,U + ,4> + \t)dWj'\t) (17c) 

where Xf , U + , V*a are probabilistic representations of position, velocity vector, and scalar 
variables, respectively. The D terms denote drift in the composition space, The B terms de- 
note diffusion, the F terms denote diffusion couplings, and the W terms denote the Wiener- 
Levy processes [29, 30]. Following Haworth and Pope [31], Dreeben and Pope [32], Colucci et 
al. [7], and Gicquel et al. [16] we consider the generalized Langevin model (GLM) and the 
linear mean square estimation (LMSE) model [26] 


dXf = U+dt + ^dW t x 
dUt = 


^- + G« (!/+-(«,» 


dxi ' '“dxkdxk 

+ + >/cwwf 

axk 


dt 




+ _ 


»* - c > (*+ - (M) + SM) 


+ 


dxkdxk 

d{<j> a ) 


dt 


dZr dW * 


(18a) 


(18b) 


(18c) 


where the variables i/i,i/ 2 , . . . are all diffusion coefficients (to be specified), and 


Gij — — U! 


\ + l Co ) 


'o 


c 

UJ — — 

k 


k 3 ! 2 1 

t — C t — k=-T(u k ,U k ) 

L\l * 


(19) 


Here u is the SGS mixing frequency, e is the SGS dissipation rate, k is the turbulent kinetic 
energy, and A l is the LES filter size. The parameters Co, C$ and C t are model constants 
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and need to be specified. The Fokker-Planck [33] equation for the diffusion process as given 
by Eq. (18) is: 

df d , r , d(p) , , d 2 (ui) ] df d [n ^ _ Ut u f] 

Tt + fe K/) = [l^ " ( " 2 “ ^ a t , |G,J M)J1 

- [-■ - + £ [c * w - <*•» 11 - 4: '* w/1 

„ ay , — a( a ,) ay . a<«_a7_ 

2 dxkdxk ^ 1/1 ^ a.c, dxidvj v 52 a*, 

1 « d(rn)d(M ay 

2 <9x,fc dxjfc 2 0 3 ^ ^Xk 0x k dvidi’a 

vs, d {<f> a ) d {<j> 0 ) d 2 / (20) 

2 dx* dxk difadrfp 

The transport equations for the filtered variables are obtained by integration of Eq. (20) 
according to Eq. (12): 


dK) _ q 

dx k 

d{ Ui ) , a («*)(«,-) _ _afr) ^ a 2 (».) giK^i) 


*•) , /fi 

dx ; V 2 


dx k dxk dxk 


(21a) 


(21b) 


^(^c) (<M 

dt dx k 




The transport equations for the second order SGS moments are 

glfe id + WlflilM = + 

dt dx k 2 dxkdxk ; 5x fc 5x fc 

+ (t-, -2 v ^ + i/ 3 ) 8ii a — - 

+ [Giir(ttt,»j) + Gjtr(tit,Uj) + C„ £ %] - g’j. " ^ 


+ (l/l - ■ 


2 dx fc dx fc v *’ 5X* 

, x d{ Ui )d{<p a ) 

'"«■ + ' /1 ^ ) + 


+ 8 W^.^) _ ^(■4I..(„ 1 _ T( „, a 

OXfc 


, , X , , M . ox dT{u k ,U t ,<p a ] 

+ [G ik T{u k , (f) a ) - C^UT[u it <p Q ) j + r(«i, b a ) ^7 


(22b) 
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dT(<j) a ,<t> 0 ) + d (Uk) 7~(^> a , <M _ tg c> 2 r(^a, <(>0) | 


dt 


dx k 

+ (v 1 - ‘IsJVxVSi + ^s 2 ) 


2 dx k dx k 

dMdM 


dx k 


dx k 


dxk dxk 


- [2 C^T^cn <f> 0 )\ + T(<f> a , Sp) + r(<^ /3 , 5„) - 


dx k 


(22c) 


A term-by-term comparison of the exact moment transport equations (Eqs. (4), (6)), with 
the modeled equations (Eqs. (21), (22)), suggests t'l = 1^2 = ^3 = ^Sj = ^s 2 = 2i/. However, 
this violates the realizability of the scalar field. A set of coefficients yielding a realizable 
stochastic model requires: V\ = 1/2 = = 2i/ and 1 / 5 , = ^s 2 = 0. That is, 

dXf = + \fTvdW* (23a) 

dUt = 


a ( P) + 2„|lM + Gij ( [/; _ (Uj » 


dan 


dxkdxk 


dt 


+ + v^dlE, 17 

d<f>t = -C > > {<f>Z - (<t>a)) dt 

The Fokker-Planck equation for this system is 


(23b) 

(23c) 
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and the corresponding equations for the moments are: 
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(25b) 

(25c) 
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Therefore, the stochastic diffusion process described by the SDEs (23) implies the following 
closure for the VSFDF: 


d_ 

dv k 


dp 

dx k 




-2v 


d 2 


dvidlpa 


dui dip a 

dx k dx k 


dx k 

v,xp\ P 


d 2 


— v- 


dui duj 
dvidvj [\dx k dx k 

d 2 T / chpa dip£ 

dx k dx k 


dip a d ip ft 


v,ipjP 

v,y) P 


d(u,)d(u,) ay i ay aw ay 

V— — — ^ r -C 0 e- - h iv- 


dx k dx k dvidvj 2 dv k dv k 


dx k dx k dvi 


■ ^ [Gi, [C ^ (V„ - (M ) /] 


(27) 


which yields the closures at the second order levels: 
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IV. NUMERICAL SOLUTION PROCEDURE 

Numerical solution of the modeled VSFDF transport equation is obtained by a Lagrangian 
Monte Carlo (MC) procedure. The basis of this procedure is the same as that used in RAS 
[34-36] and in previous FDF simulations [7, 9, 16]. In this procedure, the FDF is represented 
by an ensemble of N p statistically identical Monte Carlo particles. Each particle carries 
information pertaining to its position, X^(t), velocity, U^ n \t), and scalar value, (p^ u \t)< 
n = This information is updated via temporal integration of the SDEs. The 
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simplest way of performing this integration is via Euler-Maruyamma approximation [37]. 
For example, for Eq. (17a), 

X?(tk + i) = X t n (/ fc )+(A X (4)) n Af+(^(f fc )) n (Af) 1/2 (Cf(^)) n 

+ (F!j u (h))" (A() 1/2 (Cf(<c))“ + (F$* («)" (Al)‘ /2 (< >))" (29) 

where D x {t k ) = A(X (n) (*jt), U {n) (t k ), 4> {n) (t k )-, t k ), •••, and C(4)’s are independent stan- 
dardized Gaussian random variables. This scheme preserves the Ito character of the SDEs 

[38]. 

The computational domain is discretized on equally spaced finite difference grid points. 
These points are used for three purposes: (1) identify the regions where the statistical in- 
formation from the MC simulations are obtained, (2) perform a part ot the simulations via 
finite difference discretization coupled with the MC solver, and (3) perform a set of comple- 
mentary LES primarily by the finite difference methodology for assessing the consistency and 
convergence of the MC results. The LES procedure via the finite difference discretization is 
referred to as LES-FD and will be further discussed below. 

Statistical information is obtained by considering an ensemble of Ne computational par- 
ticles residing within an ensemble domain of characteristic length A^ centered around each 
of the finite-difference grid points. This is illustrated schematically in Fig. 1. For reliable 
statistics with minimal numerical dispersion, it is desired to minimize the size of ensemble 
domain and maximize the number of the MC particles. In this way, the ensemble statistics 
would tend to the desired filtered values: 

(a) - = — V « (n) ♦ (a) 

' /E N e “ n e - oo ' 1 

ti€.Ae A e —*Q 

r e(«» = A Y, <»>e) "(“•'>) ( ;J °) 

E n€A £ As— 0 

where a (n) denotes the information carried by n th MC particle pertaining to transport vari- 
able a. 

The LES-FD solver is based on the compact parameter finite difference scheme [39, 40]. 
This is a variant of the MacCormack scheme in which fourth-order compact differenc- 
ing schemes are used to approximate the spatial derivatives, and second-order symmetric 
predictor-corrector sequence is employed for time discretization. All of the finite difference 
operations are conducted on fixed grid points. The transfer of information from the grid 
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points to the MC particles is accomplished via a second-order interpolation. The transfer of 
information from the particles to the grid points is accomplished via ensemble averaging as 
described above. 

The LES-FD procedure determines the pressure field which is used in the MC solver. 
The LES-FD also determines the filtered velocity and scalar fields. That is, there is a 
“redundancy” in the determination of the first filtered moments as both the LES-FD and 
the MC procedures provides the solution of this field. This redundancy is actually very useful 
in monitoring the accuracy of the simulated results as shown in previous work [9, 16, 34-36]. 
To establish consistency and convergence of the MC solver, the modeled transport equations 
for the generalized second order SGS moments (Eq. (26) are also solved via LES-FD. In doing 
so, the unclosed third order correlations are taken from the MC solver. The comparison of 
all of the first and second refer moments as obtained by LES-FD with those obtained by 
the MC solver is useful to establish the accuracy of the MC solver. These simulations are 
referred to as VSFDF-C. Attributes of all of the simulation procedures are summarized in 
Table ??. In this table and hereinafter, VSFDF simulations refer to the coupled MC/LES- 
FD procedure in which the LES-FD is used for only the first order filtered variables. In 
VSFDF-C, the LES-FD procedure is used for both first and second order filtered values. 
Further discussions about the simulation methods are available in Refs. [7, 16, 34-36]. 

V. RESULTS 

A. Flows Simulated 

Simulations are conducted of a two-dimensional (2D) and a 3D temporally developing 
mixing layer involving transport of a passive scalar variable. The 2D simulations are per- 
formed to establish and demonstrate the consistency of the MC solver. The 3D simulations 
are used to assess the overall predictive capabilities of the VSFDF methodology. These 
predictions are compared with data obtained by direct numerical simulation (DNS) of the 
same layer. 

The temporal mixing layer consists of two parallel streams traveling in opposite direc- 
tions with the same speed [41-43]. In the representation below, x,y (and z) denote the 
streamwise, the cross-stream, (and the span-wise) directions (in 3D), respectively. The ve- 
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locity components along these directions are denoted by u, v , (and u>) in the x, y, (and 
z ) directions, respectively. Both the filtered streamwise velocity and the scalar fields are 
initialized with a hyperbolic tangent profiles with ( u ) = 0.5, (<f>) = 1 on the top stream and 
( u ) = —0.5, ((f)) = 0 on the bottom stream. The length L is specified such that L = 2‘ Np \ u , 
where Np is the desired number of successive vortex pairings and A u is the wavelength ol 
the most unstable mode corresponding to the mean streamwise velocity profile imposed at 
the initial time. The flow variables are normalized with respect to the half initial vorticity 

thickness, L r = M* =0 ) (£ = where ( u) L is the Reynolds averaged value of the 

filtered streamwise velocity and A 17 is the velocity difference across the layer). The reference 
velocity is U r — AU/2. 

All 2D simulations are conducted for 0 < x < 30, and —20 < y < 20. The dimensions 
are chosen to trigger the most unstable linear mode. The formation of large scale structures 
is facilitated by introducing small harmonic, phase-shifted, disturbances containing sub- 
harmonics of the most unstable mode into the stream-wise and cross-stream velocity profiles. 
For N p = 1, this results in formation of two large vortices and one subsequent pairing of 
these vortices. The 3D simulations are conducted for a cubic box, 0 < x < 60, —30 < 
y < 30, (0 < z < 60). The 3D field is parameterized in a procedure somewhat similar 
to that by Vreman et al. [44]. The formation of the large scale structures are expedited 
through eigenfunction based initial perturbations [45, 46]. This includes two-dimensional 
[42, 44, 47] and three-dimensional [42, 48] perturbations with a random phase shift between 
the 3 D modes. This results in the formation of two successive vortex pairings and strong 
three-dimensionality. 

B. Numerical Specifications 

Simulations are conducted on equally-spaced grid points with grid spacings Ax = Ay = 
Az(for 3D) = A. All 2D simulations are performed on 32 x 41 grid points. The 3D simula- 
tions are conducted on 193 3 and 33 3 points for DNS and LES, respectively. The Reynolds 
number is Re = = 50. 
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To filter the DNS data, a top-hat function of the form below is used 

3 


G(x' - x) = JJ G(x\ - x,) 


t=i 




\ x 'i ~ *i\ < 

| X \ x i | 2* 1 


( 31 ) 


No attempt is made to investigate the sensitivity of the results to the filter function [27] or 
the size of the filter [49]. 

The MC particles are initially distributed throughout the computational region. All sim- 
ulations are performed with a uniform “weight” [20] of the particles. Due to flow periodicity 
in the streamwise (and span-wise in 3D) direction(s), if the particle leaves the domain at one 
of these boundaries new particles are introduced at the other boundary with the same com- 
positional values. In the cross-stream directions, the free-slip boundary condition is satisfied 
by the mirror-reflection of the particles leaving through these boundaries. The density of 
the MC particles is determined by the average number of particles N E within the ensemble 
domain of size A E x Aj^xAe). The effects of both of these parameters are assessed to 
ensure the consistency and the statistical accuracy of the VSFDF simulations. 

All results are analyzed both “instantaneously” and “statistically.” In the former, the 
instantaneous contours (snap-shots) and scatter plots of the variables of interest are ana- 
lyzed. In the latter, the “Reynolds-averaged” statistics constructed from the instantaneous 
data are considered. These are constructed by spatial averaging over x (and r in 3D). All 
Reynolds averaged results are denoted by an overbar. 


C. Consistency and Convergence Assessments 

The objective of this section is to demonstrate the consistency of the VSFDF formulation 
and the convergence of its MC simulation procedure. For this purpose, the results via MC 
and LES-FD are compared against each other in VSFDF-C simulations. Since the accuracy 
of the FD procedure is well-established (at least for the first order filtered quantities), such 
a comparative assessment provides a good means of assessing the performance of the MC 
solution. No attempt is made to determine the appropriate values of the model constants; 
the values suggested in the literature are adopted [50] Co = 2.1, C e = 1 and C# = 1. The 
influence of these parameters are assessed in Section(V D. 
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The uniformity of the MC particle is checked by monitoring their distributions at all times, 
as the particle number density must be proportional to fluid density. The Reynolds averaged 
density field as obtained by both LES-FD and by MC are shown in Fig. 2. Close to unity 
values for the density at all times is the first measure of accuracy of simulations. Figures 
3 and 4 show the instantaneous contour plots of the filtered scalar and vorticity fields at 
several times. These figures provide a visual demonstration of the consistency of the VSFDF. 
This consistency is observed for all first order moments (filter values) without any statistical 
variability. Also, all of the first moments show very little dependence on the A E and N E 
values. So, in the presentation below we only focus on second order moments. Specifically, 
the scalar-velocity correlations are shown since all other second order SGS moments behave 
similarly. 

Figure 5-6 show the statistical variability of the results for several simulations with N E = 
40. It is observed that these moments exhibit spreads with variances decreasing as the size 
of the ensemble domain is reduced. Figures 7-10 show the sensitivity to I\ E and A E . All 
these results clearly display convergence suggested by Eq. (30). As the ensemble domain 
size decreases, the VSFDF results converge to those of LES-FD. Ideally, the LES-FD results 
should become independent of the MC results, as the latter become more reliable, i.e. when 
N E — > oo, and A E — ► 0). It is observed that best match is achieved with A^ < A/2 
and N e > 40. This conclusion is consistent with previous assessment studies on the scalar 
FDF [7, 9], and the velocity FDF [16]. All the subsequent simulations are conducted with 
A e = A/2 and N E = 40. 

D. Comparative Assessments of the VFDF 

The objective of this section is to analyze some of the characteristics of the VSFDF via 
comparative assessments against DNS data. In addition, comparisons are also made with 
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LES via the “conventional” Smagorinsky [18, 51] model: 


TL{ui,Uj) - - k 8 tl - -2 v t Sij, 


<t>) = -r t 


d (phi) 
dxi 




dx. 


dxi 


( 32 ) 
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C u = 0.04, Sc t = 1, 5 = y/SijSij and A L is the characteristic length of the filter. Thi 
model considers the anisotropic part of the SGS stress tensor a tJ = TL(ui,Uj) — 2/3 h 8 tJ . 
The isotropic components are absorbed in the pressure field. 

The predicted results via VSFDF and the Smagorinsky model are compared with those 
obtained by DNS of the same flow. For this comparison, the DNS data are filtered and are 
transposed from the original high resolution 193 3 points to the coarse 33 3 points. In the 
comparisons, we also consider the “resolved” and the “total components of the Reynolds av- 
eraged moments. The former are denoted by R(a , b) with R(a, b ) = ^(a) — (a)^ ((b) — (6)^ ; 
and the latter is r{a,b) with r(a, b) = (a - a) (b -b). n DNS, the “total” SGS components 
are directly available, while in LES they are approximated by r(a, b) « R(a, b) + r(a, b ) [44]. 
Unless indicated otherwise, the values of the model constants are (Go = 2.1, C E = 1, G^, — 1, 
but the effects of these parameters on the predicted results are assessed. 

Figure 11 show the instantaneous iso-surface of the (<f>) — 0.5 field t = 80. By this time, 
the flow has gone through several pairings and exhibits strong 3D effects. This is evident by 
the formation of large scale span-wise rollers with presence of mushroom like structures in 
streamwise planes [45]. Similar to previous results of Gicquel tt al. [16], the amount of SGS 
diffusion with the Smagorinsky model is significant, especially at initial times. Therefore, 
the predicted results are overly smooth. The Reynolds averaged values of the filtered scalar 
field at t = 80 are shown in Fig. 12 and the temporal variation of the “scalar thickness, 


MO = \vm = 0.9)1 + |»(M= 0.1)1 


(33) 


is shown in in Fig. 13. The filtered and unfiltered DNS data yield virtually indistinguishable 
results. The dissipative nature of the Smagorinsky model at initial times resulting in a 
slow growth of the layer is shown. All VSFDF predictions compare well with DNS data in 
predicting the spread of the layer. 
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Several components of the planar averaged values of the second order SGS moments are 
compared with DNS data in Figs. 14-15 for several values of the model constants. In general, 
the the VSFDF results are in better agreement with DNS data than those predicted by the 
Smagorinsky model. In this regard, therefore, the VSFDF is expected to be more effective 
than the Smagorinsky type closures for LES of reacting flows since the extent of SGS mixing 
is heavily influenced by these SGS moments [52, 53]. However, it is not possible to suggest 
“optimum” values for the model constants, except that at small C ( and values, the SGS 
energy is very large. 

Several components of the resolved second order moments are presented in Figs. 16-17. 
As expected, the performance of the Smagorinsky model is not very good as it does not 
predict the spread and the peak value accurately. The VSFDF yields reasonable predictions 
of the resolved field, except for small C<j> values. However, the total values of these moments 
are fairly independent of the model constants and yield very good agreement with DNS 
data as shown in Figs. 18-19. It is also interesting to note that while the SGS moments 
and/or the resolved moments may be over- and/or under-estimated depending on the values 
of the model coefficients, the total values of the moments are fairly independent ol these 
coefficients, at least in the range of these values as considered. But low values of C#, C' c 
are not recommended as they would result into too much SGS energy in comparison to the 
resolved energy. 

VI. SUMMARY AND CONCLUDING REMARKS 

The filtered density function (FDF) methodology has proven very effective for LES of 
turbulent reactive flows. In previous investigations, the FDF of either only the marginal 
FDF of the scalar, or that of the velocity were considered. The objective of present work is 
to develop the joint velocity-scalar FDF methodology. For this purpose, the exact transport 
equation governing the evolution of VSFDF is derived. It is shown that effects of the 
SGS convection and chemical reaction appear in a closed form. The unclosed terms are 
modeled in a fashion similar to those typically followed in PDF methods. The modeled 
VSFDF transport equation is solved numerically via a Lagrangian Monte Carlo scheme via 
consideration of a system of equivalent stochastic differential equations (SDEs). These SDEs 
are discretized via the Euler-Maruyamma approximation. The consistency of the VSFDF 
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method and the convergence of its Monte Carlo solutions are assessed in LES of a two- 
dimensional (2D) temporally developing mixing layer. This assessment is done by comparing 
the results obtained via Monte Carlo procedure with those of the finite-difference scheme 
(LES-FD) for the solution of the transport equations of the first two moments of VSFDF. 
With inclusion of the third moments from the VSFDF into the LES-FD, the consistency 
and convergence of the Monte Carlo solution is demonstrated by good agreements of the 
first two SGS moments with those obtained by LES-FD. 

The VFDF predictions are compared with those with LES results with the Smagorinsky 
[18] SGS model. All of these results are also compared with direct numerical simulation 
(DNS) data of a three-dimensional, temporally developing mixing layer. This comparison 
provides a means of examining some of the trends and overall characteristics as predicted 
by LES. It is shown that the VFDF performs well in predicting some of the phenomena 
pertaining to the SGS transport. The magnitude of the second order SGS moments as 
predicted by VSFDF is significantly larger than those predicted by the Smagorinsky model 
and are much closer to the filtered DNS results. Most of the overall flow features, including 
the mean field, the resolved and total stresses as predicted by VSFDF are in good agreement 
with DNS data. It may be possible to improve the predictive capabilities of the VFDF by two 
ways: (1) development of a dynamic procedure to determine the model coefficients, and/or 
(2) implementation of higher order closures for the generalized Langevin model parameter 
Gij (see Ref. [50]). 

Work is in progress on further development and application of VSFDF for LES of variable 
density involving exothermic chemical reactions. 
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Captions 


Table L Attributes of the computational methods. 

Figure 1. Concept of ensemble averaging. Shown are three different ensemble domains: 
1(A £ = A/2, Ne ~ 10), 2(A £ = A, N e & 40), 3(A £ = 2A, N E « 160). Black squares de- 
note the finite-difference grid points, and the circles denote the MC particles. 

Figure 2. Cross-stream variation of the Reynolds-averaged values of (p) at t=34.3: (a) 
N e = 40, (b) A E = A/2. 

Figure 3. Temporal evolution of the scalar (with superimposed vorticity iso-lines) (left) 
and the vorticity (right) fields for LES-FD, with A# = A/2 and N E = 40 at several times. 

Figure 4. Temporal evolution of the scalar (with superimposed vorticity iso-lines) (left) 
and the vorticity (right) fields for VSFDF with A E — A/2 and Ne = 40 at several times. 

Figure 5. Statistical variability of LES-FD and VSFDF simulations with N E = 40 for 
Reynolds- averaged values of r(u,^>) at t=34.4. 

Figure 6. Statistical variability of LES-FD and VSFDF simulations with N E — 40 for 
Reynolds-averaged values of at t=34.4. 

Figure 7. Cross-stream variations of the Reynolds-averaged values of r(u, <j>) (a) A# - - 
A/2, (b) A* = A, (c) As = 2A. 

Figure 8. Cross-stream variations of the Reynolds-averaged values of t(v,<p) (a) As = 
A/2, (b) As = A, (c) As = 2A. 

Figure 9. Cross-stream variations of the Reynolds-averaged values of r(u, <j>) (a) N E = 20, 
(b) N e = 40, (c) N e = 80. 

Figure 10. Cross-stream variations of the Reynolds-averaged values of r(v,<f)) (a) N E = 
20, (b) N e = 40, (c) N E = 80. 

Figure 11. Contours surface of the (<f>) field in the 3D mixing layer at t = 80. 

Figure 12. Cross-stream variations of the Reynolds averaged values of the filtered scalar 
field at t = 80. 

Figure 13. Temporal variations of the scalar thickness. 

Figure 14. Cross stream variations of some of the components of r at t = 60. 

Figure 15. Cross stream variations of some of the components of r at t = 80. 

Figure 16. Cross-stream variations of some of the components of R at t = 60. 

Figure 17. Cross-stream variations of some of the components of R at t = 80. 
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Figure 18. Cross stream variations of r at t = 60. 

Figure 19. Cross stream variations of f at t : = 80. 
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This Appendix is published as Ref. [48]. 
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Abstract. An overview is presented of recent developments and contribu- 
tions in large eddy simulation (LES) of turbulent reactive flows. The foun- 
dation of some of the recently proposed subgrid scale (SGS) closures for 
such simulations is presented, along with a discussion of their capabilities 
and limitations. The scope of the review is limited to physical modeling. 
In doing so, only issues pertaining to additional complexities caused by 
chemical reactions are discussed. That is, the challenges associated with 
“general” LES of non-reactive flows are not considered, even though all 
of these challenges are indeed present (and in most cases are a lot more 
complex) in reactive flows. It is recognized that numerical algorithms and 
computational procedures play a significant role in (any) LES. However, 
this review does not deal with these issues except for cases wherein the 
actual numerical-computational methodology is directly coupled with the 
procedure by which LES is conducted. The SGS closure based on the re- 
cently developed “filtered density function” (FDF) method is described in 
a greater detail. This is due to more familiarity of this reviewer with this 
closure; it does not imply that other closures are less effective. 


1. Introduction 

In the late 1980’s, I was preparing a review article on large scale numerical 
simulations of turbulent reactive flows. The intent was to provide a sur- 
vey of the contributions made to both direct numerical simulation (DNS) 
and large eddy simulation (LES). However, when that article was finally 
published (Givi, 1989), its content was heavily biased towards DNS. This 
was not intentional, it just reflected the state of progress on LES of turbu- 
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lent combustion at that time. But with all of the enthusiasm for DNS in 
the combustion community, the limitations of such simulations were well 
recognized (even with the most optimistic predictions of growth in super- 
computer technology). It was also clear that the future of large scale sim- 
ulations of practical turbulent reacting flows would heavily depend on the 
development of LES. Therefore, it was quite easy to predict that LES would 
receive significant attention in computational predictions of turbulent re- 
acting flows in the 1990’s and into the next (present) century. 

Now, at the time of writing this article (Summer 2001), while struggling 
to meet the deadline for its submission(!) I am not surprised by the extent 
of the contributions in developing subgrid scale (SGS) models or by the 
magnificent work on LES of a variety of turbulent reacting flow systems. 
In fact, I admit that the rate of these developments has been a lot faster 
than my capability to absorb, or in some cases even follow, the details of 
the proposed methodologies. In addition, the page-limit restrictions under 
which this is being prepared, preclude describing the details of the wide 
variety of currently available closures; similarly, citation of the relevant 
references cannot even be done exhaustively. Fortunately, many aspects of 
SGS closures and LES of reacting turbulence have recently been discussed in 
several excellent tutorial and review articles (Cook and Riley, 1998a; Candel 
et a/., 1999; Bilger, 2000; Branley and Jones, 2000; Menon, 2000; Peters, 
2000; Pope, 2000; Luo, 2001; Poinsot and Veynante, 2001). Therefore, in 
the present review I concentrate on some of the major issues related to my 
area of research within this field. 

2. Starting Equations 

Large eddy simulation involves the use of the spatial filtering operation 
(Sagaut, 2001) 


/ +°° 

<2(x',*)£(x',x)<fx', (1) 

-oo 

where Q denotes the filter function of width Aq, and (Q(x, t))e represents 
the filtered value of the transport variable Q(x,t). In variable density flows 
it is convenient to consider the Favre filtered quantity, 

(Q(x,t)) L =(pQ)i/{p)e. (2) 

We consider spatially & temporally invariant and localized filter functions, 
£7(x', x) = G(x' — x) with the properties G(x) = G{— x), and G(x)dx = 

1. Moreover, we only consider “positive” filter functions for which all the 
moments x m G(x)dx exist for m > 0. 
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To set the framework, we consider the transport equations of chemically 
reacting flows. To isolate the effects of chemical reaction in the simplest 
way, we consider single-phase (gaseous) combustion in a low Mach number 
flow with negligible radiative heat transfer and viscous dissipation. We also 
assume that Newton’s law of viscosity, Fourier’s law of heat conduction and 
Fick’s law of mass diffusion are applicable. Therefore, the primary transport 
variables are the density p, the velocity vector i = 1,2,3 along the 
X{ direction, the pressure p, the species’ mass fractions Y a , and the total 
specific enthalpy h . All of the mass fractions and the enthalpy are grouped 
into the scalar array <£(x, t) = [<f >\ , . . . 4> a \ = [Yi , F 2 , . . . , Yjv 3 , h\ of size 

a — N s + 1 where N s denotes the total number of species. Application of 
the filtering operation to the equations of continuity, momentum, enthalpy 
(energy) and species mass fraction equations gives 

d{p)t d(p) e (ui) L _ f .v 

~dT + dxi ~ °’ 


9{p)e(uj) L d_(p)e{u,) L (uj)L 
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dt 


dx{ 
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dxj dxi dxi 
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d{p)l{<t>a)L d{p)t{Uj)L{<j>a)L _ 

dt dx i dx{ 


dM f 
dxi 


+ (pS a )ij 


( 5 ) 


where t represents time, and the filtered reaction source terms are denoted 
by {pSa)i = {p)i{S a )i. The viscous stress tensor and the mass/heat fluxes 
are denoted by r 4 y, and Jf* , respectively. At low Mach numbers and heat 
release rates, by neglecting the viscous dissipation and thermal radiation 
the source terms in the enthalpy equation can be assumed to be negligible. 
Thus, S a = S a (<f>). The terms T %3 - (p)t((uiu j)l - ( Ui)i{uj) L ) and Aff = 
(p)e({u>i<f>a)L ~ ( u i ) l(4>q)l) denote the SGS stress and the SGS mass flux, 
respectively. Equations (3)-(5) are coupled through the equation of state. 


3. Closure Methodologies 

For non-reacting flows the SGS closure is associated with Tij and M-* 
(Canuto, 1994; Ciofalo, 1994; Lesieur and Metais, 1996). In reacting flows, 
an additional model is required for the filtered reaction rate This 

modeling is the subject of primary concern in this review. 

One of the first contributions in LES of reactive flows, similar to that in 
LES of non-reacting flows, was made in atmospheric sciences (Schumann, 
1989). In this work, the effects of SGS scalar fluctuations (as appear in the 
chemical source term) are assumed negligible, i.e. (S a (<f>))L ~ 5'a({^)L)- 
This assumption is compatible with that made in some of the more recent 
contributions (Boris et a/., 1992; Fureby and Grinstein, 1999), in which it 



84 


PEYMAN GIVI 


is argued that all of the essential SGS contributions are included in the 
numerical discretization procedure. 

Modeling of the scalar fluctuations has been the subject of broad in- 
vestigations in Reynolds averaged simulations (RAS) for over five decades, 
resulting in a variety of closure strategies (Libby and Williams, 1980; Libby 
and Williams, 1994). Within the past 10 years or so, almost all of these 
closures have been considered for LES. Examples: the eddy-break up mod- 
els (Fureby and Lofstrom, 1994; Candel et a/., 1999), moment methods 
(Frankel et a/., 1993), the flamelet concept (Cook et a/., 1997; Cook and 
Riley, 1998b; De Bruyn Kops et al., 1998; DesJardin and Frankel, 1998; 
DesJardin and Frankel, 1999; Pitsch and Steiner, 2000; Ladeinde et a/., 
2001), the linear eddy model (LEM) (McMurtry et a/., 1992; Menon and 
Calhoon, 1996; Kim et al . , 1999; Menon, 2000), the conditional moment 
method (CMM) (Bushe and Steiner, 1999; Steiner and Bushe, 2001), and 
many others (Sykes et a/., 1992; Galperin and Orszag, 1993; Smith and 
Menon, 1996; Im et al . , 1997; McGrattan et a/., 1998; Thibaut and Candel, 
1998; Battaglia et a/., 2000; Collin et a/., 2000). In addition, several of the 
closures previously developed for LES of non-reacting flows, have been ex- 
tended for use in reacting flow simulations (DesJardin and Frankel, 1998; 
Jaberi and James, 1998). 

The probability density function (PDF) methods have proven particu- 
larly useful in RAS (O’Brien, 1980; Pope, 1985; Dopazo, 1994; Fox, 1996; 
Pope, 2000). The systematic approach for determining the PDF is by means 
of solving its transport equation. An alternative approach is based on 
assumed methods in which the shape of the PDF is specified a priori. 
This has been pursued in several studies in most of which it is assumed 
that the thermo-chemical variables depend only on the mixture fraction, 
e.g. infinitely fast reaction, equilibrium chemistry. Therefore, the PDF is 
univariate (Madnia and Givi, 1993; Cook and Riley, 1994; Reveillon and 
Vervisch, 1996; Branley and Jones, 1997; Jimenez et al., 1997; Mathey and 
Chollet, 1997; DesJardin and Frankel, 1998; DesJardin and Frankel, 1999; 
Forkel and Janicka, 2000; Kempf et al., 2000). For LES of non-equilibrium 
reactive flows, it is necessary to assume the joint PDF of multi-scalars 
(Frankel et al., 1993). Consistent with popular methods of generating uni- 
variate (Leemis, 1986) and multivariate (Johnson and Kotz, 1972) distribu- 
tions, all of the assumed SGS scalar PDFs in the contributions cited above 
are based on the first and the second order moments. The PDFs generated 
in this way offer sufficient flexibility and are affordable for large scale sim- 
ulations. However, it is now well understood that the “true” PDF strongly 
depends on the actual physics of mixing in a given flow condition (Jaberi 
et al., 1996). Therefore, there is a need to determine such PDFs in a more 
systematic manner. 
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The “filtered density function” (FDF) methodology introduced by Pope 
(1990) provides the framework for fundamental developments of the PDF 
based SGS closures. This method provides a means of determining the PDF 
from its own transport equation. For the scalars’ array 0(x, £) the FDF, 
denoted by Pl, is defined as (Pope, 1990) 

/ +oo 

C [Vs </>(*', 0] G( x ' -x)dx, (6) 

-OO 


a 

< [V>, <£(X, 0] = II &[$ a ~ <M X , *)], (7) 

a=l 

where <5 denotes the delta function and 0 denotes the composition domain 
of the scalar array. The term ([</> — 0(x, t)] is the “fine-grained” density 
(Lundgren, 1967; O’Brien, 1980; Pope, 1985; Dopazo, 1994). In variable 
density flows, it is convenient to consider the “filtered mass density func- 
tion” (FMDF), denoted by Fl, as 

/ +oo 

p(x', t)C [V>, ^(x', <)] G(x' - x)dx!. (8) 

•OO 


The integral property of the FDF and FMDF is such that 

/ +oo y+o o 

= 1, / F L (if>;x,t)diJ> = (p(x,t)) e . (9) 

-oo J — OO 

For further discussions, it is useful to define the mass weighted conditional 
filtered mean of the variable Q(x,/), 


f-£ p(x', t)Q(x', tX [V>, <£(x', 0] G{x! - x)dx' 


{Q(x,t)\il>)e = 


J — OO 


F L (ij>;x,t) 


( 10 ) 


Therefore, when Q can be completely described by the compositional vari- 
able, i.e. Q(x,t ) = Q(<j>(x,t)), we have ( Q(x,t)\i>)t = Also, 



(Q(x,t)\tl>) e F L (rp\x,t)d^ = {p(x,t)) e {Q{x,t))L- 


( 11 ) 


The transport equation for F/,(-0 ; x, t) is obtained by multiplying the trans- 
port equation for the fine grained density by the filter function G(x' — x) 
and integrating over x' space (Gao and O’Brien, 1993; Colucci et al., 1998; 
Reveillon and Vervisch, 1998; Jaberi et al., 1999; Jaberi, 1999; Zhou and 
Pereira, 2000; Tong, 2001), 
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dfjAjh x, t) d[{uj(x,t) |- 0 )^i ? L(' 0 ;x,O] 

dt dx{ 


d[S Q {Tl>)F L (xlr,x,t)\ 

di> a 


+ 


dlp a 


1 dJf 


p(4>) dx 


L IV’\ F L (ip;x, 

1 f l 


<) 


( 12 ) 


The first term on the RHS is due to chemical reaction and is in a closed 
form. This demonstrates the primary advantage of the FDF methodology. 
However, the SGS convection (the second term on the LHS) and SGS mix- 
ing (the second term on the RHS) must be modeled. One of the most chal- 
lenging issues in FDF is associated with closure of the mixing term. This 
has been the subject of broad investigations in PDF modeling (Pope, 1985; 
Pope, 2000). In Eq. (12) the effects of mixing are displayed through the 
“conditional expected diffusion” of the scalars, but can also be represented 
in the form of the “conditional expected dissipation” (O’Brien, 1980; Pope, 
1985). The closure for this can be via any of the ones currently in use in 
PDF methods (Pope, 2000). In the absence of a clearly superior model, the 
linear mean square estimation (LMSE) model (O’Brien, 1980) has been 
used in almost all of previous LES based on FDF (Colucci et a/., 1998; 
Jaberi et a/., 1999; Garrick et a/., 1999; James and Jaberi, 2000; Zhou and 
Pereira, 2000). With Jf — — this model is 


— U- 1 -— 

dtpa l\ pdxi 



= d \ d{F L l(p) t ) ] 

dxi l dx{ J 

r\ 

+ [Hm(0a - (0a)L)Ct,] ,(13) 


where fl m (x,2) is the “frequency” of mixing within the subgrid and must 
be modeled. The convective term can be modeled as 

{ui\ip)tF L = (wj)l-Cl - 

where 7 * is the SGS diffusion coefficient and must be specified. Equation 
(14) is in accord with that often used in conventional LES (Moin et a/., 
1991; Canuto, 1994; Ciofalo, 1994; Lesieur and Metais, 1996). With this 
formulation, obviously the resolved hydrodynamic field must be determined 
by other means. This problem can be circumvented by considering the joint 
velocity- scalar FMDF, 



/ +00 

p(x', t)i [v, u(x', t), -0, 0(x', I)] G(x > - x)dx', (15) 

-OO 
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3 a 

i [v, u(x, t), -0, <£(x, 2)] = n 6[v k - U k (x , <)] n 0]» ( 16 ) 

k=l cx=l 

where v denotes the composition domain of the random velocity vector, and 
f[v, u(x, £), *0, <£(x, t)] is the fine-grained velocity-scalar density. Most re- 
cent work in this regard consider the transport of the velocity FDF (VFDF) 
(Gicquel, 2001) and the joint velocity-scalar FDF (VSFDF) (Drozda, 2001). 
The operational procedure is similar to that developed previously for PDF 
methods (Pope, 1985; Pope, 1994; Pope, 2000). 

The closure problems as noted above are not particular to the FDF; 
all of the other schemes require similar modelings. For example, in the 
limit of equilibrium chemistry all of the statistics of the reacting fields 
are related to those of the mixture fraction. The FMDF of the mixture 
fraction can be obtained from the solution of Eq. (12) with S = 0. So, 
there is still a need for modeling of the mixing term. Even in cases where 
the FMDF is assumed, its distribution is parameterized with the low order 
moments of the mixture fraction. As indicated above, the first two moments 
are typically used for this parameterization. Therefore, there is a need for 
closure of the “total SGS dissipation” as appears in the second moment 
(SGS variance) equation. Several means of dealing with this closure problem 
are available (Girimaji and Zhou, 1996; Pierce and Moin, 1998; Cook and 
Bushe, 1999; Jimenez et a/., 2001; De Bruyn Kops and Riley, 2001). 

The above problem is a bit more complex when the SGS chemical reac- 
tion is assumed to be in the “flamelet” regime (Peters, 2000). In this case, 
even with the one-dimensional flamelet model, the thermo-chemical vari- 
ables are parameterized by the mixture fraction and its rate of dissipation 
(Cook et a/., 1997; Cook and Riley, 1998b; De Bruyn Kops et a/., 1998; 
DesJardin and Frankel, 1998; Cook and Riley, 1998b). Therefore, there is a 
need for a priori specification of the joint FDF of the mixture fraction and 
its dissipation. A review of different methods of dealing with this issue is 
available (Cook and Riley, 1998a). Equation (12) with S’ = 0 indicates that 
there is a dependence between the FDF of the mixture fractions and the 
conditional expected diffusion (and the conditional expected dissipation). 
This dependency is not considered in most previous contributions, but is 
the subject of current investigations (DesJardin et a/., 2001). 

Modeling of the conditional expected dissipation is also required in the 
conditional moment method (Bushe and Steiner, 1999; Steiner and Bushe, 
2001). This issue has been recognized at the early stages of developments of 
CMM in RAS (Bilger, 2000). With this model, the conditional filtered mean 
values of the thermo- chemical variables (LHS of Eq. (10)) are obtained 
by their modeled transport equation. This is obviously computationally 
less demanding that solving the FDF transport equation. But in order to 
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determine the actual filtered quantities, the distribution of the mixture 
fraction FDF must be specified. 

An important issue in regard to FDF is associated with the numeri- 
cal solution of its transport equation. The Lagrangian Monte Carlo scheme 
(Pope, 1985) has proven particularly useful for this purpose. In this scheme, 
the FDF is represented via an ensemble of computational elements or par- 
ticles. Transport of these particles and the change in their properties are 
modeled by a set of stochastic differential equations (SDEs) (Soong, 1973). 
The diffusion process (Gardiner, 1990) has proven effective for this pur- 
pose. The coefficients in the Langevin equation governing this process are 
set in such a way that the resulting Fokker-Planck equation (Risken, 1989) 
is equivalent to the FDF transport equation. Therefore, the Monte Carlo 
solution of the SDEs represent the solution of the FDF in the probabilistic 
sense. This procedure has proven successful for simulating PDF in a variety 
of systems (Grigoriu, 1995). However, one must be careful in performing 
stochastic simulations in conjunction with modern CFD solvers. Many of 
the advanced discretization routines developed for solving deterministic dif- 
ferential equations may not be applicable, or may have to be significantly 
modified to be suitable for solving SDEs (Kloeden and Platen, 1995). 

Implementation of LEM is also based on stochastic representation of 
the flow. In its original development in RAS (Kerstein, 1988), the processes 
of molecular diffusion, chemical reaction and turbulent convection are con- 
sidered separately. This is achieved by a reduced one-dimensional (linear) 
description of the scalar field, which makes it possible to resolve the flow 
scales even for flows with relatively high Reynolds, Schmidt and Damkohler 
numbers. The interpretation of the one-dimensional domain is dependent 
on the particular flow under consideration. In this way, the processes of 
molecular diffusion and chemical reaction are taken into account exactly, 
but the effects of convection are modeled. This is achieved by “random 
rearrangement” (or stirring) events in such a way that the displacements 
of fluid elements result in a diffusivity equal to the “turbulent diffusivity.” 
For LES, this procedure is followed within each of the computational cells, 
and stirring is performed to yield the desired SGS diffusivity. Menon and 
colleagues have made extensive use of LEM for LES of a wide variety of 
reacting flows. A recent review is available (Menon, 2000). 
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